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Introduction 

The  purpose  of  our  work  is  to  develop  and  test  a  new  kind  of  imaging  method 
we  call  "palpation  imaging”.  We  expect  that  palpation  imaging  will  be  a 
useful  tool  for  improving  the  discrimination  between  benign  and  malignant 
breast  tumors.  The  scope  of  the  effort  in  our  first  two  years  of  funding  was  to 
implement  a  newly-developed  algorithm  for  motion  tracking  on  a  commercial 
ultrasound  imaging  system  and  to  begin  testing  that  new  imaging  system. 
The  algorithm  provides  images  of  the  mechanical  strain  induced  in  tissue  by 
pressing  the  ultrasound  transducer  against  the  skin  surface.  These  images 
are  produced  at  substantially  real-time  frames  rates  with  normal  ultrasound 
B-mode  and  strain  images  displayed  side-by-side.  The  algorithm  is  fully 
integrated  into  the  commercial  system  and  requires  no  system  modifications. 
Breast  exams  performed  on  volunteers  have  shown  that  the  palpation  imaging 
study  is  almost  identical  to  the  standard  clinical  breast  ultrasound  exam. 
With  the  patients  examined  to  date,  we  find  distinct  strain  image  sequences 
for  cysts,  fibroadenomas,  and  invasive  ductal  carcinomas.  A  comparison  of 
lesion  area  measured  in  B-mode  versus  strain  images  appears  to  provide  good 
discrimination  between  benign  and  malignant  lesions  with  the  data  obtained 
so  far.  This  work  is  summarized  in  a  manuscript  published  in  Ultrasound  in 
Medicine  and  Biology  29(3);  427-35,  2003  (attached  as  Appendix  1  in  the 
annual  report  last  year). 

The  plans  for  the  third  and  final  year  of  the  proposed  effort  included 
a  larger,  more  well-defined  study  of  the  performance  of  palpation  imaging 
for  breast  lesion  classification.  We  planned  to  scan  enough  human  subjects 
to  obtain  at  least  100  patients  with  biopsy  of  surgery-proven  lesion  t3q)e. 
Unfortunately,  there  has  been  a  significant  shortfall  in  the  human  subject 
recruitment  for  this  project.  This  is  due  to  a  lack  of  cooperation  of  the 
cUnician  co-investigators  involved  in  this  project.  The  shortage  of  medical 
staff  and  an  increased  demand  to  provide  clinical  care  has  left  these  co¬ 
investigators  with  limited  time  for  research  activities.  Those  changes  as  well 
as  an  outstanding  opportunity  resulted  in  the  Principle  Investigator,  Timothy 
Hall,  leaving  the  University  of  Kansas  Medical  Center  (KUMC)  for  a  position 
as  a  tenured  Professor  in  the  Medical  Physics  Department  at  the  University 
of  Wisconsin-Madison  (UW).  The  benefits  of  this  change  of  environment 
are  wide-ranging  and  include  a  large  department  of  investigators  pursuing 
advances  in  medical  imaging,  a  vast  pool  of  graduate  students  and  post¬ 
doctoral  candidates  to  work  with,  and  a  medical  staff  actively  involved  in 


our  research. 

Expenditure  of  funds  at  KUMC  ceased  18.Jan.20035  two  days  after  Pro¬ 
fessor  Hall  started  work  at  UW.  Although  work  has  continued  on  this  project, 
at  both  KUMC  and  UW,  no  additional  funds  have  been  spent  against  the 
grant.  We  have  been  working  toward  transferring  this  grant  to  UW  since  that 
time.  It  appears  that  the  process  is  nearly  complete.  Associated  with  the 
grant  transfer  we  have  requested  a  one-year  no-cost  extension  of  this  grant 
period.  This  will  provide  ample  time  to  complete  that  proposed  study. 


Body 

Four  of  the  six  proposed  Tasks  (Tasks  1-3  and  6)  have  been  addressed  so  far 
within  this  project.  The  overall  effort  has  been  very  successful  and  has  drawn 
attention  of  the  elasticity  imaging  research  community  around  the  world.  We 
are  extremely  excited  with  our  results  to  date  and  are  very  encouraged  with 
the  potential  for  this  new  kind  of  imaging  system.  Below  is  a  description  of 
the  approved  Tasks  and  the  progress  toward  achieving  each  goal. 

Task  1.  Implement  real-time  palpation  imaging  on  a  commercial  sonog¬ 
raphy  system  (months  1-4)*  o,)  Program  the  imaging  system  digital  signal 
processors  to  estimate  strain  from  consecutive  image  frames,  b)  Develop  a 
user  interface  for  controlling  the  data  acquisition  and  processing  parameters, 
c)  Test  the  system  using  existing  phantoms  and  laboratory  fixtures  with  mo¬ 
torized  motion  control  to  determine  the  penalty  for  using  fixed-point  versus 
floating  point  calculations. 

This  task  has  been  completed,  as  described  in  the  year-2  report.  We  are 
grateful  to  SMSUG  for  their  technical  assistance  in  this  effort.  The  soft¬ 
ware  we  have  developed  is  fully  functional,  and  has  gone  through  multiple 
improvements.  The  current  version  displays  3  cm  wide  by  3.5  cm  deep  im¬ 
age  pairs  at  about  seven  frames  per  second  (73.5  cm^/sec).  That  frame  rate 
is  sufficiently  fast  to  provide  feedback  to  the  hand-eye  coordination  system 
to  allow  manipulating  the  conditions  of  tissue  compression  to  consistently 
obtain  high-quality  strain  images.  The  manuscript  that  describes  the  imple¬ 
mentation  of  this  algorithm  and  initial  testing  for  displacement  variance  is 
published  (Ultrasonic  Imaging  24(3):  161-176,  2002). 

Task  2.  Develop  data  acquisition  techniques  that  provide  high-quality 
palpation  images  without  the  use  of  fixtures  (months  3-12):  a)  Implement 
techniques  that  mimic  those  used  for  phantom  imaging  including  small  hand- 
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held  fixtures  to  restrict  motion  perpendicular  to  the  image  plane,  b)  Test  those 
techniques  in  phantoms  and  compare  target  contrast-to-noise  ratio  for  labora¬ 
tory  (large  motorized  fixtures j  controlled  motion)  versus  clinical  (small  hand¬ 
held  fixtures  ^  restricted  motion)  systems,  c)  Test  the  clinical  systems  (real¬ 
time  palpation  imaging  with  small  fixtures)  using  anthropomorphic  breast 
phantoms,  d)  Test  the  clinical  systems  on  volunteer  patients  with  palpable 
breast  abnormalities,  e)  Modify  the  data  acquisition  techniques  to  eliminate 
the  need  for  fixtures  to  restrict  motion  while  maintaining  image  quality,  f) 
Measure  conspicuity  of  breast  lesions  to  assess  the  relative  merit  of  different 
data  acquisition  techniques. 

This  task  has  been  completed,  as  described  in  the  year-2  report.  Iter¬ 
ating  between  what  we  learned  with  freehand  scanning  of  phantoms,  then 
in  vivo  breasts,  and  back  to  phantoms,  we  found  that  the  key  to  obtaining 
sequences  of  high-quality  strain  images  (high  contrast-to-noise  and  high  sim¬ 
ilarity  from  frame  to  frame)  with  freehand  scanning  is  high  frame  rate.  The 
scanning  technique  for  breast  palpation  imaging  is  nearly  identical  to  that 
used  in  typical  breast  sonography  with  compression.  The  subject  hes  on  her 
back  with  her  arm  behind  her  head  and  the  ultrasound  transducer  is  pressed, 
by  hand,  toward  the  chest  wall.  The  real-time  image  feedback  allow  manip¬ 
ulation  of  the  compression  conditions  while  scanning.  This  manipulation  is 
essential  for  obtaining  long  sequences  of  low  noise  strain  image  data  (often 
100  image  frames). 

Task  3.  Implement  high-quality  palpation  imaging  algorithm  on  a  com¬ 
mercial  sonography  system  and  perform  preliminary  tests  of  image  quality 
(months  7-17):  a)  Program  the  commercial  sonography  system  to  calculate 
high-resolution,  low-noise  palpation  images  as  quickly  as  possible  over  a  large 
region  of  interest,  b)  Develop  a  user  interface  that  allows  manipulation  of 
the  image  formation  algorithm  for  the  trade-off  between  spatial  resolution 
and  image  noise,  c)  Test  the  high-quality  algorithm  on  (geometrically)  sim¬ 
ple  and  anthropomorphic  phantoms  using  the  modified  data  acquisition  tech¬ 
niques  (Task  2.e).  d)  Use  the  real-time  palpation  imaging  technique  to  locate 
the  desired  region  of  interest  and  obtain  sonographic  data  with  the  appro¬ 
priate  pre-  and  post-compression  for  volunteer  patients  with  palpable  breast 
abnormalities,  and  Task  6.  Investigate  the  use  of  novel  techniques,  such  as 
harmonic  imaging  and  spatial  quadrature,  for  improved  information  content 
in  palpation  images  (months  18-36). 

This  area  of  work  has  been  our  dominant  effort  since  the  first  year  of 
support.  As  described  in  the  year-1  report,  we  can  estimate  displacement 
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and  strain  accurately  in  inhomogeneous  media  (tissue)  with  greater  success 
than  originally  anticipated.  The  key  to  this  success  is  real-time  imaging  of 
strain  that  guides  manipulation  of  the  conditions  of  compression. 

As  described  in  previous  reports,  we  have  implemented  three  displace¬ 
ment  and  strain  estimation  algorithms  on  the  Elegra.  The  simplest  of  these 
is  the  algorithm  that  runs  in  real-time  for  data  acquisition.  This  uses  our 
patented  adaptive  search  strategy  to  predict  the  deformation  based  on  previ¬ 
ous  displacement  estimates.  This  adaptive  search  reduces  the  computational 
load  for  strain  estimation  by  more  than  a  factor  of  100.  The  algorithm  also 
decimates  the  data,  thereby  reducing  the  number  of  displacement  estimates 
in  a  given  region  of  interest.  This  algorithm  is  used  primarily  for  data  ac¬ 
quisition  during  freehand  scanning  and  is  a  compromise  between  frame  rate 
and  image  quality.  To  halt  data  acquisition,  the  system  ‘freeze’  button  is 
pressed,  just  as  in  normal  sonography,  and  the  echo  data  are  available  for 
on-line  post-processing.  This  same  (real-time)  strain  imaging  algorithm  can 
be  used  to  reprocess  that  data  as  the  user  scrolls  through  image  memory,  or 
one  of  the  other  algorithms  can  be  used.  The  second  algorithm  is  identical 
to  the  first  except  that  the  data  is  not  decimated.  Estimating  displacement 
with  higher  data  density  reduces  the  displacement  variance  in  the  adaptive 
search  strategy  resulting  in  a  slight  improvement  in  image  quality.  The  third 
algorithm  operates  on  the  full  data  field  (no  decimation)  and  does  not  use  the 
adaptive  search  strategy,  instead  using  a  full  2-D  search.  This  approach  is 
much  more  computationally  intensive,  requiring  about  one  second  per  strain 
image  frame,  but  produces  images  that  lack  the  displacement  error  accmnu- 
lation  that  can  result  from  the  adaptive  search  strategy. 

We  are  continuing  to  investigate  options  to  further  improve  image  qual¬ 
ity,  as  described  in  the  year-2  annual  report.  These  advances  are  expected 
to  significantly  improve  on-line  post-processed  images.  As  a  result  of  that 
effort  it  has  become  clear  that  subjective  visual  assessment  of  image  ‘qual¬ 
ity’  and  lesion  conspicuity  are  not  sufficient  to  quantify  differences  among 
motion  tracking  algorithms.  As  a  result,  we  are  working  on  a  quantitative 
method  for  measuring  the  accuracy  of  motion  tracking  in  tissues  in  addition 
to  the  estimates  of  displacement  estimate  error  variance  that  are  commonly 
reported.  This  work  is  very  early  in  development,  and  is  not  expected  to  be 
completed  during  this  project  period.  But,  the  approach  shows  great  promise 
for  comparing  the  relative  performance  of  various  motion  tracking  algorithms 
and  parameter  selection  for  a  specific  algorithm. 

To  date  we  have  scanned  56  volunteers  at  KUMC.  Some  of  these  patients 
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had  multiple  lesions,  and  some  were  scanned  on  repeat  occasions  (perpendic¬ 
ular  image  planes,  repeat  visits,  different  transducers,  etc).  As  reported  last 
year,  the  excitement  from  the  potential  for  palpation  imaging  and  the  de¬ 
lay  at  KUMC  in  obtaining  a  sufficient  number  of  patients  to  adequately  test 
this  technology  resulted  in  Siemens  contracting  with  two  other  institutions  to 
test  our  technology  (Charing  Cross  Hospital  in  London,  UK  and  Mayo  Clinic, 
Rochester,  MN).  Those  sites  have  been  quite  productive  in  their  trials.  To 
date  the  Charing  Cross  group  has  scanned  over  200  breast  patients  and  the 
Mayo  group  has  scanned  over  100  patients.  Each  found  that  they  became 
proficient  in  their  scanning  technique  within  10-15  patients,  and  they  find 
that  palpation  imaging  adds,  on  average,  about  five  additional  minutes  to 
the  normal  breast  ultrasound  examination.  The  key  to  strain  imaging,  as 
in  standard  sonography,  is  training.  With  real-time  feedback  of  the  current 
strain  images  being  acquired,  the  key  is  interpreting  the  data  and  making  cor¬ 
rections  to  the  technique,  if  necessary.  In  summary,  the  tools  are  adequate, 
the  key  to  success  is  training. 

As  described  in  the  previous  annual  reports,  we  have  demonstrated  ev¬ 
idence  for  changing  strain  image  contrast  in  fibroadenomas  as  the  surface 
pressure  increases.  Support  for  that  observation  continues  to  build,  and  that 
behavior  continues  to  appear  unique  to  fibroadenomas.  The  previous  an¬ 
nual  reports  also  described  the  frame-to-frame  variability  of  strain  images 
for  cysts.  Although  those  observations  are  still  true  of  the  cysts,  our  empha¬ 
sis  in  data  acquisition  has  been  on  solid  lesions.  We  also  provided  evidence 
in  previous  reports  describing  the  frame-to-frame  similarity  in  strain  images 
of  invasive  ductal  carcinomas  and  evidence  demonstrating  that  the  apparent 
size  of  a  carcinoma  in  strain  images  is  larger  than  that  seen  in  the  B-mode 
images.  Results  reported  in  year-1  were  for  measurements  performed  by  the 
PI.  A  more  extensive  study,  described  below,  was  reported  at  the  Radiological 
Society  of  North  America  meeting  in  November  2002. 

Garra,  et  al,,  (Radiology  202:  79-86,  1997)  suggested  that  the  width  of 
a  carcinoma  in  a  strain  image  is  typically  larger  than  that  measured  in  a 
B-mode  image.  Our  results  support  that  observation,  and  apparently  extend 
its  diagnostic  utility,  as  described  in  the  previous  annual  reports.  To  begin 
designing  the  experiments  for  Tasks  4  and  5,  we  asked  five  observers  (the  PI 
and  four  ultrasound  clinicians)  to  view  each  sequence  of  side-by-side  B-mode 
and  strain  images  from  the  KUMC  data.  Prom  each  sequence,  each  observer 
had  to  choose  an  image  pair  that  represented  the  typical  view  of  the  lesion 
in  both  B-mode  and  strain  images.  That  image  pair  was  then  displayed 
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and  the  observer  had  to  trace  the  lesion  boundary  in  the  B-mode  image  and 
measure  the  lesion  height  and  width.  The  image  pair  was  then  re-displayed 
and  the  observer  traced  the  lesion  in  the  strain  image  and  measure  the  lesion 
height  and  width.  The  combined  average  results  for  lesion  measurements 
were  shown  in  the  year-2  report.  A  receiver  operating  charateristic  (ROC) 
analysis  of  those  measurements  resulted  in  an  overly-optimistic  diagnostic 
accuracy  of  0.983±0.028— nearly  perfect  performance.  Although  this  is  a 
very  encouraging  results,  it  cannot  be  trusted  for  several  reasons  including 
the  small  size  of  the  patient  pool,  correlations  among  the  data  (multiple 
lesions  in  a  single  patient  included  in  the  study),  and  the  limited  variety  of 
lesion  types  (invasive  ductal  carcinomas,  fibroadenomas  and  cysts). 

However,  this  work  demonstrates  that  we  are  following  a  solid  research 
plan.  By  increasing  the  size  of  the  study  population,  and  comparing  the 
change  in  diagnostic  accuracy  when  using  standard  sonography  versus  sonog¬ 
raphy  plus  palpation  imaging,  we  can  test  the  utility  of  palpation  imaging  as 
a  diagnostic  tool  with  ROC  analysis  and  multi-observer  measurements. 

Our  current  plans  are  unchanged  from  those  in  the  proposal  and  current 
protocol,  except  for  the  change  of  institution  and  the  use  of  data  from  outside 
this  funded  project.  We  are  well  prepared  to  complete  the  final  two  tasks 
of  the  study:  Tasks  4-  Acquire  sonograms,  elastograms,  and  mammograms 
on  patients  with  suspicious  lesions  that  are  either  palpable  or  detectable  with 
sonography  or  mammography  (months  18-32;  >100  patients) and 
Task  5.  Compare  diagnosis  of  breast  lesions  with  and  without  the  use  of 
palpation  imaging,  (months  30-36;  using  images  obtained  in  Task  4)- 

Key  Research  Accomplishments 

•  The  motion  tracking  algorithm  has  been  implemented  on  the  Siemens 
SONOLINE  Elegra  and  displays  B-mode  and  strain  images  side-by-side 
at  about  seven  frames  per  second. 

•  We  found  that  the  key  to  obtaining  high-quality  in  vivo  strain  images 
of  the  breast  is  to  form  the  images  in  “real-time,”  that  is,  fast  enough  to 
provide  the  hand-eye  coordination  system  sufficient  feedback  to  control 
the  conditions  of  tissue  deformation. 

•  We  have  developed  scanning  techniques  for  acquiring  strain  image  data 
from  in  vivo  breasts  that  provides  reproducible  results.  The  technique 
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is  nearly  identical  to  typical  breast  sonography  with  compression,  so  it 
is  relatively  easy  to  learn  and  clinicians  are  likely  to  adopt  the  practice. 

•  There  is  no  significant  difference  between  the  displacement  variance 
(strain  image  noise)  for  freehand  scanning  versus  motorized  compres¬ 
sion  for  strain  imaging. 

•  We  have  found  that  the  frame-to-frame  strain  patterns  fi-om  various 
breast  abnormalities  appears  to  be  unique  to  the  abnormahty.  For  ex¬ 
ample,  the  fluid  within  cysts  appears  to  be  ‘stirred’  when  deformed  in 
palpation  iihaging,  and  that  ‘stirring’  causes  the  rf  echo  signal,  and 
therefore  the  strain  image,  to  decorrelate  rapidly  with  time.  We  have 
also  found  that  the  strain  contrast  for  fibroadenomas  is  not  constant 
with  compression;  at  very  low  surface  pressure  fibroadenomas  are  stiffer 
than  their  surroundings  and  provide  high  negative  contrast.  With  in¬ 
creased  surface  pressure  that  contrast  is  often  significantly  reduced; 
Invasive  ductal  carcinomas  maintain  a  high  negative  contrast  at  all 
surface  pressures  tested. 

•  We  have  found  that  by  comparing  the  area  of  a  lesion  measured  on 
the  standard  ultrasound  B-mode  image  with  area  measured  on  a  strain 
image,  benign  lesions  have  nearly  equal  area  on  both  modalities  but 
invasive  ductal  carcinomas  are  significantly  larger  on  the  strain  image 
than  in  B-mode. 

•  The  diagnostic  criterion  of  comparing  lesion  size  in  B-mode  versus 
strain  images  appears  to  be  a  strong  criterion  among  a  group  of  five 
observers.  It  must  be  noted  that  this  was  a  relatively  srnall  data  set 
with  some  lesions  that  were  easily  classified  without  the  benefit  of  pal¬ 
pation  ima.ging  results.  This  criterion  will  be  further  tested  in  Tasks  4 
and  5. 

Reportable  Outcomes 
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Conclusions 

i 

Our  initial  research  plan  included  the  development  and  initial  testing  of  a 
method  for  real-time  imaging  of  mechanical  strain  in  tissue  and  is  proceeding 
as  planned.  Our  success  in  this  effort  far  exceeds  our  anticipated  results.  The 
phantom  studies  proposed  proved  far  too  simple  to  evaluate  the  the  merit 
of  this  technique,  but  the  in  vivo  results  demonstrate  the  reproducibility  of 
strain  imaging.  Measurements  of  lesion  size  in  strain  images  compared  to  B- 
mode  images  appears  to  be  a  sensitive  diagnostic  criterion  for  discriminating 
malignant  from  benign  lesions.  Further  improvements  in  the  motion  tracking 
algorithm  will  make  strain  imaging  more  robust  and  increase  the  confidence 
of  the  sonographer.  Our  results  to  date  indicate  that  real-time  palpation 
imaging  has  the  potential  to  significantly  improve  the  diagnosis  of  breast 
abnormalities.  This  new  tool  runs  as  a  software  application  on  an  existing 
clinical  sonography  system  and  is  therefore  easily  distributed  broadly  when  it 
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is  an  appropriate  product.  Our  software  is  computationally  intensive,  but  it 
is  expected  that  future  sonography  systems  from  all  manufacturers  will  have 
sufficient  computational  capacity  to  support  this  application.  Wide-spread 
use  would  likely  follow,  if  our  current  results  are  representative.  We  look 
forward  to  continued  success. 
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IN  VIVO  REAL-TIME  FREEHAND  PALPATION  IMAGING 
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Abstract — ^Previous  experience  with  laboratory  fixtures  and  ofiT-line  processing  of  elasticity  data  showed  that 
problems  occurring  in  data  acquisition  often  resulted  in  poor  elasticity  image  quality.  A  system  for  real-time 
estimation  and  display  of  tissue  elastic  properties  using  a  clinical  ultrasonic  imaging  system  has  been  developed. 
A  brief  description  of  that  system  and  the  initial  clinical  tests  of  that  system  are  reported.  Experience  with 
real-time  freehand  elasticity  imaging  shows  that  images  with  high  contrast-to-noise  ratios  are  consistently 
obtained.  Images  of  breast  lesions  were  acquired  with  freehand  palpation  using  standard  linear-array  ultrasound 
(US)  transducers.  Results  in  volunteer  patients  show  that  high-quality  elasticity  images  are  easily  obtained  from 
in  vivo  breast  studies.  The  key  element  to  successful  scanning  is  real-time  visual  feedback  of  B-mode  and  strain 
images  that  guide  the  patient  positioning  and  compression  direction.  Results  show  that  individual  images  of  axial 
strain  in  tissues  can  be  quite  misleading,  and  that  a  ^movie  loop’’  of  side-by-side  B-mode  and  strain  images 
provides  significantly  more  information.  Our  preliminary  data  suggest  that  the  strain  image  sequences  for 
various  breast  patholo^es  are  unique.  For  example,  strain  images  of  fibroadenomas  lose  contrast  with  increasing 
precompression,  but  those  of  invasive  ductal  carcinoma  have  high  negative  contrast  (dark  relative  to  "normal” 
tissue)  for  a  wide  range  of  precompression.  In  addition,  a  comparison  of  the  lesion  area  measured  in  B-mode  vs. 
strain  images,  for  a  representative  image  from  the  sequence,  appears  to  be  a  sensitive  criterion  for  separating 
invasive  ductal  carcinoma  from  cyst  and  fibroadenoma.  (E-mail:  thall@wisc.edu)  €>  2003  World  Federation  for 
Ultrasound  in  Medicine  &  Biology. 

Key  Words:  Ultrasound,  Tissue  characterization.  Elasticity,  Palpation,  Elastography. 


INTRODUCTION 

The  potential  for  improving  the  qualitative  nature  of 
palpation  by  imaging  quantitative  measures  of  tissue 
viscoelasticity  has  generated  a  great  deal  of  interest 
world-wide.  Our  initial  efforts  focused  on  modeling  dis¬ 
placement  and  strain,  developing  algorithms  for  dis¬ 
placement  and  strain  estimation,  and  testing  those  tech¬ 
niques  in  phantoms  and  in  kidneys  in  vitro  (see,  for 
example,  Chaturvedi  et  al.  1998a,  1998b;  Hall  et  al. 
1997;  Insana  et  al.  1997;  Zhu  et  al.  1998).  Significant 
effort  was  expended  on  developing  high-order  motion 
estimators  for  tracking  fine-scale  motion.  However,  little 
data  were  available  to  investigate  the  need  or  utility  of 
the  high-order  motion-estimation  techniques  for  in  vivo 
imaging  of  tissues. 

The  first  report  testing  the  utility  of  strain  imaging  in 
breast  lesion  imaging  ((jarra  et  al.  1997)  clearly  demon¬ 
strated  that  strain  imaging  had  merit  in  breast  lesion 
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discrimination.  The  data-acquisition  system  employed  a 
modified  mammography  compression  paddle  and,  there¬ 
fore,  was  limited  in  the  lesion  locations  that  could  be 
studied.  Also,  only  (nonreal-time)  static  strain  images 
were  available.  In  that  report,  Garra  et  al.  (1997)  de¬ 
scribed  a  set  of  criteria  applied  to  evaluate  strain  imaging 
combined  with  normal  B-mode  imaging.  Among  those 
criteria  were  lesion  visibility,  relative  brightness,  lesion 
margin  regularity,  lesion  margin  definition,  lesion  size 
(lateral  and  axial)  and  B-mode  image  measurements  rel¬ 
ative  to  strain  image  and  pathology  measurements. 
Among  their  findings,  they  noted  that  benign  lesions 
have  about  the  same  width  on  B-mode  and  strain  images, 
but  that  the  height  measurement  had  lower  confidence 
due  to  axial  blurring  in  strain  image  formation  and  dif¬ 
ficulty  in  determining  lesion  boundaries  with  shadows 
due  to  high  attenuation.  Fibroadenomas  typically  had 
heterogeneous  stiffness;  cancers  were  uniformly  stiffer 
than  their  surroundings  in  all  but  one  case. 

The  purpose  of  the  present  study  was  to  test  the 
utility  of  performing  strain  imaging  in  real-time  on  a 
commercial  ultrasound  (US)  imaging  system  and  to  test 
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one  of  the  strain  image  criteria  described  by  Garra  et  al. 
(1997)  with  data  acquired  from  this  new  system.  Results 
demonstrate  the  value  in  real-time  side-by-side  display 
of  B-mode  and  strain  images  for  guiding  data  acquisition 
and  data  interpretation.  Comparisons  among  various  le¬ 
sion  types  studied  in  vivo  show  a  significant  difference  in 
the  strain  image  sequence  for  fibroadenomas,  cysts  and 
carcinoma,  and  help  to  explain  some  of  the  difficulties  in 
data  interpretation  described  by  Garra  et  al.  (1997).  Our 
results  are  generally  consistent  with  those  found  by 
Garra  and  colleagues,  but  the  differences  we  found  in 
carcinoma  size  in  B-mode  and  strain  images  are  greater, 
and  aH  lesions  found  in  sonography  or  mammography, 
whether  palpable  or  not,  were  visible  with  our  tech¬ 
niques.  These  results  will  help  to  guide  future  strain¬ 
imaging  data  acquisition  and  provide  further  evidence 
for  the  potential  of  elasticity  imaging  in  breast  lesion 
discrimination. 

MATERIALS  AND  METHODS 

The  motion-tracking  algorithm,  its  implementation 
on  the  commercial  clinical  US,  imaging  system  and 
performance  measurements  in  experiments  with  phan¬ 
toms  are  reported  elsewhere  (Zhu  and  Hall  2002).  The 
essential  details  are  included  here  for  the  convenience  of 
the  reader. 

Strain  image  formation 

A  2-D  block-matching  algorithm,  based  on  the  sum- 
squared  difference  (SSD)  method,  was  used  for  motion 
tracking  in  our  implementation.  With  this  method,  mo¬ 
tion  is  tracked  by  searching  for  a  kernel  of  data  from  the 
precompression  radio  fi^uency  (RF)  echo  data  in  a  2-D 
search  region  of  the  postcompression  RF  echo  field,  A 
fixed  kernel  size  (five  A-lines  by  11  RF  samples)  was 
used  with  both  the  7.5L40  (with  and  without  tissue 
harmonic  imaging,  THI)  and  the  VFX13-5  linear  arrays 
for  the  system  employed  (Siemens  SONOLINE  Elegra, 
Issaquah,  WA),  The  data  are  temporally  sampled  at  36 
MHz  with  a  lateral  beam  spacing  of  200  /im.  Therefore, 
the  kernel  size  corresponds  to  about  about  one  half  the 
area  of  the  2-D  pulse-echo  US  point  spread  function  for 
the  7.5L40  array  pulsed  at  7.2  MHz.  A  small  kernel  was 
chosen  because  the  assumption  of  rigid  body  motion  is 
increasingly  accurate  as  the  size  of  the  kernel  is  de¬ 
creased,  and  because  spatial  resolution  is  expected  to 
improve  with  smaller  kernels.  No  attempt  to  optimize  the 
kernel  size  was  pursued  in  this  study.  Kernel  size  opti¬ 
mization  will  likely  be  task-dependent  and  will  be  ad¬ 
dressed  in  future  work. 

Data  were  processed  on  the  image  processor  sub¬ 
system  of  the  Elegra.  This  subsystem  hosts  two  Texas 
Instruments  TMS320C80  processors.  To  reduce  the  com¬ 
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putational  load  (required  to  achieve  real-time  frame 
rates),  an  adaptive  search  strategy  was  developed  that 
reduces  the  size  of  the  required  search  region  in  perform¬ 
ing  the  SSD  block  matching.  Displacements  are  esti¬ 
mated  (for  real-time  imaging)  with  kernels  that  are  sep¬ 
arated  by  16  RF  samples  center-to-center  (no  spatial 
overlap  of  displacement  estimates).  Displacements  are 
estimated  row-by-row,  and  the  prior  row  of  estimates  are 
used  to  predict  displacements  in  the  current  row,  allow¬ 
ing  the  search  region  to  be  reduced  to  within  one  RF 
sample  in  each  direction  of  the  predicted  displacements. 
The  use  of  predicted  displacements  results  in  correlated 
displacement  errors,  and  an  error  detection  and  correc¬ 
tion  scheme  was  also  implemented.  Strain  is  estimated 
from  the  displacement  data  using  a  linear  regression 
technique  similar  to  that  described  by  Kallel  and  Ophir 
(1997).  For  real-time  imaging,  linear  regression  is  per¬ 
formed  with  a  24-sample  window  (about  4  mm)  that  is 
incremented  one  displacement  sample  for  each  strain 
estimate.  The  resulting  algorithm  displays  streaming  B- 
mode  and  strain  images  side-by-side  at  about  seven 
frames  per  second  and  stores  the  full  sequence  of  I-Q 
(analytic  form  of  the  RF)  echo  data  at  full  system  bus 
speed.  The  stored  data,  which  were  acquired  at  a  higher 
frame  rate  than  the  real-time  display,  can  then  be  online 
postprocessed  with  the  same  displacement  algorithm  or 
other  algorithms  (not  reported  here),  the  size  and  location 
of  the  subregion-of-interest  (SROI)  can  be  adjusted,  the 
grey-scale  mapping  can  be  modified,  etc.,  and  the  results 
displayed  frame-by-frame  or  as  a  cine  loop. 

Although  the  additional  degrees  of  freedom  of  mo¬ 
tion  allowed  with  freehand  compression,  compared  with 
motorized  compression,  were  expected  to  result  in  an 
increase  in  displacement  estimate  error  variance,  that 
variance  is  only  slightly  higher  for  freehand  compression 
(Zhu  and  Hall  2002).  Given  that  the  contrast  and  reso¬ 
lution  of  the  strain-imaging  system  do  not  depend  on  the 
method  of  tissue  deformation,  displacement  estimate  er¬ 
ror  variance  (resulting  in  strain  image  noise)  is  the  dom¬ 
inant  image  quality  parameter  that  will  differ  with  the 
two  methods  of  tissue  deformation.  With  equal  applied 
strain,  motorized  and  freehand  compression  have  com¬ 
parable  strain  image  noise.  However,  the  frame-to-frame 
strain  is  not  constant  with  fi^ehand  compression  (as 
described  below).  This  results  in  fiame-to-frame  vari¬ 
ability  in  strain  image  quality  with  freehand  scanning. 

Small  (eg,,  2,4-mm  diameter)  isoechoic  spherical 
targets  in  a  phantom  are  considerably  easier  to  locate  and 
scan  freehand  than  with  motorized  compression.  The  size 
of  spherical  targets  measured  from  the  resulting  strain 
images  is  very  close  to  their  true  dimensions  (both  height 
and  width,  see  Zhu  and  Hall  2002),  so  both  linear  and 
area  measurements  in  strain  images  in  vivo  should  be 
accurate. 
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The  grey-scale  mapping  of  these  strain  images  con¬ 
forms  with  the  de  facto  standard  of  mapping  pixels, 
representing  small  strains  dark  and  large  strains  bright. 
Typically,  a  sinusoidal  (freehand)  compress/release  cy¬ 
clic  deformation  was  used  and  the  acquired  data  con¬ 
tained  one  or  more  complete  cycles.  At  the  top  and 
bottom  of  the  stroke,  there  is  little  motion  and,  occasion¬ 
ally,  there  were  brief  hesitations  during  the  motion.  A 
consequence  of  freehand  scanning  is  that  the  fiame-to- 
frame  strain  is  not  constant.  To  compensate  for  that 
variability  in  average  strain,  the  grey-scale  for  individual 
strain  images  (in  the  sequence  of  B-mode  and  strain 
images)  was  automatically  adjusted  to  minimize  bright¬ 
ness  flicker.  With  this  scaling,  strain  images  that  had  an 
average  strain  of  less  than  0.15%  were  set  to  black.  These 
flames  rarely  occur,  but  are  most  common  at  the  top  and 
bottom  of  the  cyclic  compression.  Frames  with  greater 
than  0.15%  compression  were  encoded  as  follows: 


where  sfk,  0  is  the  encoded  strain  value  at  position  {k,  /), 
and  and  are  the  maximum  and  minimum  strain 
values  in  the  frame,  and  roundQ  is  a  function  that  roxmds 
to  the  nearest  integer.  After  scaling,  the  results  are  en¬ 
coded  into  an  8-bit  display  range.  This  simple  scaling 
provides  a  reasonably  constant  strain  image  brightness 
through  the  compression  cycle,  in  the  absence  of  signif¬ 
icant  displacement  estimation  errors.  Large  displacement 
estimation  errors  sometimes  occur  and  this  automatic 
grey-scale  mapping  can  be  dominated  by  erroneous  dis¬ 
placement  estimates. 


Patient  scanning 

All  patients  provided  informed  consent  consistent 
with  the  protocol  approved  by  the  Human  Subjects  Com¬ 
mittee  (Institution^  Review  Board)  at  The  University  of 
Kansas  Medical  Center.  Patient  scans  were  performed  in 
a  manner  consistent  with  a  normal  breast  US  examina¬ 
tion;  the  breast  was  scanned  with  the  patient  (typically) 
in  the  supine  position  with  her  ipsilateral  arm  behind  her 
head.  When  the  breast  lesion  was  located,  SROI  was 
chosen  that  would  avoid  inappropriate  data  (lungs,  areas 
of  lost  transducer  contact,  etc.)  and  the  transducer  was 
pressed  toward  the  chest  wall  at  a  steady  rate  in  an  effort 
to  achieve  about  1-1.5%  compression  fi:ame-to-frame. 
Subregion  selection  also  typically  excluded  the  retro¬ 
mammary  fat  layer  and  the  chest  wall.  The  soft  fat,  the 
stiff  muscle  and  the  slipping  boundary  between  these 
layers  can  also  dominate  the  dynamic  grey-scale  map¬ 
ping.  In  some  cases,  for  example,  when  scanning  lateral 


lesions  in  large  (D-cup)  breasts,  the  patient  was  rolled 
slightly  to  her  contralateral  side  so  Aat  gravity  would 
flatten  the  breast  tissue  in  the  region  to  be  scanned.  A 
small  plate  (approximately  45-mm  wide,  90-mm  long) 
was  sometimes  attached  to  the  transducer  body  to  extend 
the  compression  surface  in  an  effort  to  provide  a  more 
uniform  stress  field  and  to  control  motion  perpendicular 
to  the  image  plane.  The  compress/release  cycle  was 
repeated  for  relatively  large  (>  10%)  compression  range, 
while  watching  the  B-mode  image.  The  compression 
motion  was  adjusted  by  changing  the  compression  direc¬ 
tion  or  patient  position  until  there  was  nearly  uniaxial 
compression  with  minimal  elevation  motion.  With  this 
achieved,  the  strain-imaging  software  was  enabled  to 
evaluate  the  quality  of  the  sequence  of  strain  images.  If 
a  large  sequence  (^  30  frames)  of  strain  images  had 
good  image  quality  (relatively  high  contrast-to-noise  ra¬ 
tio)  and  high  frame-to-frame  similarity,  the  data  acqui¬ 
sition  was  frozen,  the  image  sequence  stored,  and  the 
cine  feature  of  this  software  was  used  to  review,  post¬ 
process  and  select  images  to  record.  If  the  compression 
rate  was  too  slow,  resulting  in  low  frame-average  strain, 
the  interframe  skip  was  adjusted  to  increase  the  strain 
between  frame  pairs  used  in  displacement  and  strain 
estimation,  as  suggested  by  Lubinski  et  al.  (1999).  Rep¬ 
resentative  results  obtained  when  scanning  a  3-mm  cyst 
are  shown  in  Fig.  1.  The  average  strain  per  frame  (Fig. 
lb)  suggests  nearly  ideal  compression  rate  in  this  case. 
Consecutive  frames  were  paired  for  displacement  esti¬ 
mation  when  analyzing  this  sequence  of  data.  The  cu¬ 
mulative  strain  in  the  sequence  (Fig.  Ic),  obtained  by 
summing  the  strain  in  consecutive  finmes,  demonstrates 
that  about  an  18%  compression  range  was  achieved  in 
this  study.  The  initial  value  on  the  ordinate  axis  in  each 
plot  was  set  to  zero  and  (the  frame-average)  strain  was 
accumulated  (Fig.  Ic)  from  that  starting  point.  If,  for 
example,  there  was  a  net  compression  of  the  tissue  in  the 
first  frame  pair,  the  initial  average  strain  was  negative 
and  the  first  value  on  the  cumulative  frame  plot  was 
negative.  Data  acquisition  was  frozen  by  the  operator 
when  an  acceptable  image  sequence  was  acquired.  The 
starting  and  ending  points  of  that  sequence  could  be  at 
any  point  in  the  cyclic  motion.  Therefore,  the  precom¬ 
pression  at  the  starting  point  is  a  random  value  between 
the  minimum  and  maximum  precompression. 

RESULTS 

One  of  the  most  promising  uses  of  straiii  imaging  is 
differentiation  among  breast  lesions.  To  date,  we  have 
successfully  scanned  39  patients  with  over  175  patient 
scans  (multiple  lesions  per  patient,  repeat  visits,  perpen¬ 
dicular  scan  planes,  different  transducers  and  THI).  Data 
from  only  29  of  these  patients  are  included  in  this  study 
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Fig.  1.  Data  obtained  by  freehand  scanning  of  a  breast  cyst  in 
vivo,  (a)  A  B>mode  and  strain  image  pair  obtained  for  frame  5 1 
in  the  sequence.  The  white  box  in  the  B-mode  image  defines  the 
ROI  for  strain  imaging,  (b)  The  average  strain  in  the  ROI;  (c) 
the  cumulative  strain  in  the  ROI.  The  arrows  in  (b)  and  (c) 
indicate  that  a  sequence  of  frames  was  acquired  with  nearly 
equal  average  strain  in  each  frame  but  with  varying  cumulative 
strain  (precompression). 


(data  from  patients  with  surgical  scars  or  closely  spaced 
lesions  were  excluded,  as  described  below).  Among 
these  29  patients,  19  fibroadenomas,  29  cysts  and  7 
carcinomas  were  included  in  this  study.  Ten  of  the  fi¬ 
broadenomas,  at  least  five  of  the  cysts  and  four  of  the 
carcinomas  were  palpable.  All  fibroadenomas  were  ei¬ 
ther  pathologically  proven  or  had  been  stable  under  ra¬ 
diological  investigation  for  more  than  1  year.  All  carci¬ 
nomas  were  pathologically  proven.  Both  simple  and 
complex  (hemorrhagic)  cysts  were  included.  Fibroade¬ 
nomas,  cysts  and  invasive  ductal  carcinomas  have  dis¬ 
tinctive  behavior  in  their  strain  image  under  cyclic  com¬ 
pression,  as  detailed  below. 

One  of  the  key  tests  was  to  show  that  strain  images 
are  reproducible,  both  within  an  image  sequence  and  on 
repeat  acquisitions.  The  question  of  reproducibility 
within  an  image  sequence  is  addressed  for  each  lesion 
type  (fibroadenoma,  cyst  and  invasive  ductal  carcinoma) 
below.  Figure  2  shows  results  of  repeating  the  freehand 
in  vivo  elasticity  imaging  on  the  same  patient,  A  skilled 
sonographer  can  acquire  a  sufficiently  similar  B-mode 
image  of  an  ROI  with  repeat  scanning.  However,  obtain¬ 
ing  a  similar  strain  image  requires  that  the  ROI  be  found, 
and  the  direction  of  compression/release  relative  to  the 
lesion  and  chest  wall  be  the  same  in  the  two  studies.  No 


m 


Fig.  2.  Two  image  pairs  from  the  same  patient  on  repeat  visits. 
The  patient  has  a  palpable  fibroadenoma  that  measures  about 
16  mm  X  11  mm.  (a)  The  strain  image  acquired  with  the 
VFX13-5  array  during  the  first  visit,  which  is  very  similar  to  (b) 
that  acquired  with  the  7.5L40  array  2  weeks  later. 


special  effort  was  used  to  obtain  similar  images,  but  Fig. 
2  demonstrates  that  the  strain  patterns  in  these  images  are 
very  similar. 

Too  few  independent  samples  of  each  lesion  type 
have  been  observed  to  make  strong  statistical  statements 
regarding  each  criterion  described  by  Gana  et  al.  (1997). 
The  following  descriptions  state  our  qualitative  observa¬ 
tions  to  date,  in  the  hope  of  guiding  the  scanning  and 
image-evaluation  techniques  used  by  others  in  future 
studies.  In  particular,  the  automatic  scaling  of  grey-scale 
values  precludes  quantitative  statements  of  (and  homo¬ 
geneity  of)  relative  stiffness. 

A  total  of  37  B-mode  and  strain  image  sequences 
were  acquired  from  19  unique  fibroadenomas  among  9 
patients.  One  of  our  most  significant  findings  is  that  there 
was  an  obvious  (subjective,  visual)  loss  in  strain  image 
contrast  for  14  of  these  19  fibroadenomas  (27  of  37 
image  sequences,  6  of  9  patients,  average  age  44  years). 
Fibroadenomas  typically  have  negative  contrast  (are 
more  stiff  than  their  surroundings)  at  low  precompres¬ 
sion  and  lose  contrast  (stiflfhess  becomes  more  like  their 
surroundings)  as  they  are  compressed.  An  example  of 
this  is  illustrated  in  Fig.  3.  The  largest  negative  contrast 
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Fig.  3.  Data  obtained  by  freehand  scanning  of  a  fibroadenoma 
in  vivo,  (a)  The  average  strain  in  the  ROI  for  each  frame 
suggests  a  slow  compression  rate.  The  interframe  skip  was 
increased  to  pair  every  fourth  frame  in  analyzing  this  sequence 
of  data,  (b)  The  cumulative  strain  in  the  sequence  shows  that  a 
20%  compression  range  was  achieved.  Images  of  fibroadeno¬ 
mas  acquired  with  equivalent  precompression  [equivalent  cu¬ 
mulative  strain,  frames  (c)  75  and  (e)  90,  (d)  52  and  (f)  106]  are 
similar,  but,  because  of  the  nonlinear  stress-strain  relationship 
of  tissue,  strain  images  acquired  at  different  precompression 
can  have  significantly  different  contrast. 


occurs  with  minimal  precompression.  Figure  3  c  and  e 
shows  that  the  transducer  is  barely  in  contact  with  the 
skin  surface  at  diis  precompression.  Comparing  these 
with  Fig.  3d  and  f  demonstrates  that  images  of  fibroad¬ 
enomas  acquired  within  a  sequence  at  equivalent  pre¬ 
compression  are  very  similar.  However,  as  precompres¬ 
sion  changes,  strain  contrast  changes  (comparing  Fig.  3c 
or  e  with  d  or  f. 

Both  simple  and  complex  cysts  were  included  in  this 
study,  to  investigate  the  possibility  that  fluid-filled  cysts, 
regardless  of  “echogenicity,”  have  a  common  behavior  in 
cyclic  strain  images.  The  hope  was  that  strain  images 
might  help  to  differentiate  complex  cysts  from  solid 
lesions  that  lack  evidence  of  blood  flow. 

The  frame-to-frame  variability  of  strain  images  of 
cysts  is  more  complicated  than  that  observed  with  fibro¬ 
adenomas.  A  total  of  39  B-mode  and  strain  image  se¬ 
quences  were  acquired  from  29  unique  cysts  among  15 
patients.  A  very  soft  bottom  layer  in  (he  interior  of  the 
cyst  was  observed  in  7  of  29  cysts  (11  of  12  image 
sequences  of  those  cysts).  That  layer  might  be  due  to  a 
sediment  inside  the  cystic  fluid.  Repositioning  the  patient 


Fig.  4.  Data  obtained  by  freehand  scanning  of  two  breast  cysts 
in  vivo,  (a)  The  cumulative  strain  in  the  sequence  shows  that  a 
20%  compression  range  was  achieved.  Strain  images  of  the 
interior  of  cysts,  unlike  those  of  fibroadenomas,  are  not  neces¬ 
sarily  similar  when  acquired  with  similar  precompression  (cu¬ 
mulative  strain),  (b),  (c)  and  (d)  Images  from  frames  8, 25  and 
92,  respectively,  were  acquired  with  similar  precompression. 
Although  strain  images  can  vary  smoothly  from  frame  to  frame, 
decorrelation  of  the  signals  within  the  cysts  results  in  strain 
images  that  vary  significantly  over  the  compression  cycle. 


might  have  allowed  us  to  confirm  that  conjecture,  but 
that  test  was  not  performed.  The  interior  echoes  within 
the  cysts  rapidly  decorrelate  with  compression.  As  a 
result,  the  apparent  strain  in  the  lesions  varies  with 
compression,  but  (hat  compression-dependent  strain  im¬ 
age  contrast  is  very  different  from  that  observed  for 
fibroadenomas.  If  the  incremental  average  strain  from 
one  strain  image  to  the  next  is  small  (<  0,5%),  the  strain 
image  brightness  (on  a  pixel-by-pixel  basis)  changes 
gradually,  regardless  of  precompression.  If  the  incremen¬ 
tal  average  strain  is  not  small  (S:  1%),  then  the  local 
brightness  within  the  cyst  varies  rapidly  and  (seemingly) 
unpredictably.  The  typical  behavior  of  the  strain  pattern 
in  cysts  is  demonstrated  (as  well  as  can  be  with  static 
images)  in  Fig.  4.  Unlike  the  behavior  observed  with 
fibroadenomas,  frames  with  equivalent  precompression 
might  have  very  different  apparent  strain  within  the  cyst 
Overall,  a  cyst  can  be  either  relatively  stiff,  as  if  it  were 
a  distended  balloon,  or  relatively  soft. 

A  total  of  21  B-mode  and  strain  image  sequences 
were  acquired  from  7  unique  carcinomas  among  6  pa¬ 
tients.  All  carcinomas  studied  so  far  were  invasive  ductal 
carcinomas  and  all  but  one  were  highly  suspicious  of 
carcinoma,  based  on  mammogram  and  sonogram  results. 
This  is  by  far  the  most  commonly  diagnosed  breast 
cancer.  Relatively  small  lesions  (^  2  cm)  have  high 
negative  contrast  (stiff)  in  a  background  of  “normal” 
breast  tissue,  regardless  of  precompression.  An  example 
of  this  is  shown  in  Fig.  5.  The  exception  to  this  occurs  for 
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(b)  (c) 


Fig.  5.  Data  obtained  by  freehand  scanning  of  invasive  ductal 
carcinoma  in  vivo,  (a)  The  cumulative  strain  in  the  sequence 
shows  that  a  8%  compression  range  was  achieved,  (b)  Images 
in  frame  4  acquired  at  low  precompression,  are  very  similar  to 
(c)  those  in  frame  18  acquired  at  much  higher  precompression. 
Some  differences  in  the  strain  in  “normal”  surrounding  tissue 
are  seen,  but  vaiy  smoothly  from  frame  to  frame. 


(«) 


Fig.  6.  B-mode  and  strain  images  of  frbroadenomas  with  their 
perimeter  traced  in  the  B-mode  image;  that  tracing  also  in  the 
strain  image  for  comparison.  A  fibroadenoma  measuring  (a) 
71.4  mm^  in  B-mode  and  75.0  mm^  in  strain,  (b)  88,7  mm^  in 
B-mode  and  102  mm^  in  strain,  (c)  27.2  mm^  and  21.5  mm^  in 
B-mode  and  26.8  mm^  and  21.7  mm^  in  strain,  respectively. 
The  B-mode  tracing  is  a  reasonably  good  approximation  to  that 
on  the  strain  image. 


very  large  lesions  where  little,  if  any,  healthy  tissue  is 
included  in  the  strain  image. 

One  of  the  criteria  that  Garra  et  al.  (1997)  foimd  to 
be  most  useful  in  differentiating  between  benign  and 
malignant  lesions  was  the  relative  size  of  the  lesion  in 
B-mode  vs.  strain  images.  To  compare  lesion  size  in  the 
two  imaging  modalities,  we  transferred  die  I-Q  echo  data 
to  an  off-line  computer  for  further  analysis.  The  strain 
images  were  reprocessed  using  the  same  displacement 
estimation  algorithm  as  that  implemented  on  the  Elegra. 
Off-line  processing  used  a  16-sample  (<  3  mm)  window 
for  strain  estimation  instead  of  the  24-sample  window 
used  on  the  Elegra  (higher  axial  resolution).  Movie  loops 
of  the  side-by-side  B-mode  and  strain  image  pairs  (avi 
files)  were  created  to  view  the  motion  of  the  lesion  in  the 
B-mode  image  and  the  resulting  strain  image.  A  repre¬ 
sentative  frame  was  selected  that  showed  the  “typical” 
strain  image  for  that  lesion  and  the  B-mode  image  was 
displayed,  allowing  the  lesion  boundary  to  be  traced.  The 
boundary  in  the  B-mode  image  excluded  the  capsule  of 
the  lesion.  The  lesion  width  (and  height)  were  estimated, 
based  on  the  traced  lesion  perimeter,  as  the  maximum 
dimension  perpendicular  (and  parallel)  to  the  acoustic 
beam.  The  tracing  and  measurement  process  was  then 
repeated  with  the  strain  image  from  that  same  frame.  The 
boundary  traced  in  the  strain  image  was  the  location  of 
the  steepest  (visual)  gradient  in  strain.  High  negative- 
contrast  images  were  chosen  for  fibroadenomas.  All  trac¬ 
ings  were  performed  by  the  first  author,  and  most  bound¬ 
aries  were  very  easily  identified.  In  some  cases,  for 
example,  in  the  atypical  fibroadenoma  shown  in  Fig.  3, 


an  experienced  clinician  assisted  in  tracing  the  boundary. 
Example  images  for  a  fibroadenoma,  cysts  and  an  inva¬ 
sive  ductal  carcinoma  are  shown  in  Figs.  6,  7  and  8. 

It  is  intriguing  to  examine  the  relative  size  of  these 
lesions,  comparing  their  width,  height  and  area  as  mea¬ 
sured  in  B-mode  and  strain  images.  Figure  9a  and  b 
shows  plots  of  the  width  and  height  of  these  three  lesion 
types  as  measured  in  B-mode  and  strain  images.  Figure 
9c  shows  plots  of  a  similar  comparison  of  the  total  area 


(4  ft) 


(c) 

Fig.  7.  B-mode  and  strain  images  of  cysts  with  their  perimeter 
traced  in  the  B-mode  image;  that  tracing  also  in  tihe  strain 
image  for  comparison.  A  cyst  measuring  (a)  139  mm^  in 
B-mode  and  145  mm^  in  strain,  (b)  102  mm^  in  B-mode  and 
84.1  mm^  in  strain,  (c)  30.3  mm^  in  B-mode  and  32,5  mm^  in 
strain.  The  B-mode  tracing  is  a  good  approximation  to  that  on 
the  strain  image,  if  the  soft  region  at  the  bottom  of  the  cyst,  (a) 
and  (b),  were  included. 
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Fig,  8.  B-mode  and  strain  images  of  invasive  ductal  carcinomas 
with  their  perimeter  traced  in  the  B-mode  image;  that  tracing 
also  in  the  strain  image  for  comparison.  An  invasive  ductal 
carcinoma  measuring  (a)  96.1  mm^  in  B-mode  and  170  mm^  in 
strain,  (b)  22.8  mm^  in  B-mode  and  319  mm^  in  strain,  (c)  48.7 
mm^  in  B-mode  and  1 70  mm^  in  strain,  (d)  465  mm^  in  B-mode 
and  at  least  768  mm^  in  strain.  The  B-mode  tracing  is  not 
representative  of  what  would  likely  have  been  drawn  on  any  of 
these  strain  images. 


of  the  lesion  in  the  two  imaging  modes.  Table  1  shows 
that  the  width  and  height  of  benign  lesions  tend  to  be 
about  the  same  size  in  B-mode  and  strain  images  and 
carcinomas  are  larger  in  strain  images  than  in  B-mode,  as 


observed  by  Garra  et  al.  (1997),  but  the  separation  be¬ 
tween  benign  lesion  and  carcinoma  is  larger  when  we  use 
the  lesion  area. 

In  each  of  these  examples,  the  study  was  performed 
on  an  isolated  lesion.  Although  some  data  sets  contained 
more  than  one  lesion,  those  lesions  were  separated  by  at 
least  one  diameter  of  the  largest  lesion  in  the  image. 
Measurements  of  individual  lesions  in  clusters  of  lesions, 
most  frequently  observed  in  clusters  of  breast  cysts, 
proved  problematic  in  obtaining  high-quality  strain  im¬ 
age  sequences  and  in  interpreting  the  motion.  An  exam¬ 
ple  of  this,  shown  in  Fig.  10,  shows  that  it  might  be  more 
reasonable  to  study  the  cluster  as  a  group  instead  of  as 
individual  lesions.  When  lesions  are  closely  spaced,  par¬ 
ticularly  when  they  share  a  common  boundary,  the  mo¬ 
tion  due  to  compression  can  be  quite  complex,  as  ob¬ 
served  in  the  B-mode  image  sequence,  and  the  block¬ 
matching  algorithm  fails  to  track  motion  adequately.  The 
block-matching  algorithm  assumes  rigid  body  motion 
and  does  not  accurately  track  significant  rotation  or  shear 
motion.  Further,  our  current  system  acquires  data  (effec¬ 
tively)  in  a  plane,  and  significant  motion  perpendicular  to 
the  image  plane  is  lost.  A  3-D  acquisition  system  would 
be  required  to  track  significant  elevation  motion.  A  high¬ 
er-order  algorithm,  such  as  the  deformable  mesh  (Zhu  et 
al.,  1998),  would  be  required  to  track  rotation  and  shear 
motion  accurately.  At  this  stage  of  strain  image  process¬ 
ing  and  interpretation,  it  is  likely  best  to  restrict  Ae  study 
to  individual  isolated  lesions. 


(a) 


(b) 


(c) 


(d) 


Fig.  9.  Plots  comparing  the  size  of  a  lesion  traced  in  the 
B-mode  image  vs.  the  same  lesion  traced  in  a  representative 
strain  image  for  (O)  cysts,  (□)  fibroadenomas  and  (x)  invasive 
ductal  carcinomas,  (a)  The  width  and  (b)  height  comparisons, 
(c)  The  areas  for  lesions  less  than  20-mm  wide  (in  B-mode)  are 
compared,  and  (d)  the  area  ratio.  The  dashed  line  in  (a)-(c) 
represents  equal  size  measurement  in  both  images. 


DISCUSSION 

Real-time  display  of  side-by-side  B-mode  and  strain 
images  is  essential  for  guiding  the  manipulation  of 
boimdary  conditions  for  the  mechanics  experiment  that  is 
strain  imaging.  Real-time  feedback  to  the  hand-eye  co¬ 
ordination  system  allows  the  sonographer  to  manipulate 
the  compression  direction,  force  and  rate  to  obtain  high- 
quality  sequences  of  strain  images.  The  system  employs 
standard  linear-array  transducers  and  requires  no  addi¬ 
tional  fixtures  or  remote  data-acquisition  or  signal-pro- 
cessing  hardware.  It  is  fully  integrated  into  the  Elegra 
system.  A  small  plate  9  cm  X  4.5  cm  is  sometimes 
attached  to  the  face  of  the  transducer  to  extend  file 
compression  surface.  This  was  most  useful  when  scan¬ 
ning  benign  lesions  that  tended  to  move  in  elevation 
when  compressed. 

Our  results  show  significantly  different  strain-image 
sequences  for  each  lesion  type  studied.  Although  the 
three  lesion  types  reported  here  do  not  include  all  those 
found  in  breasts,  they  represent  the  most  common  clin¬ 
ically  observed  breast  lesions.  It  was  found  that,  to  ap¬ 
preciate  the  differences  among  lesion  types,  and  to  de¬ 
termine  the  “typical”  strain  image  for  a  given  lesion,  it 
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Table  1 .  Results  of  measurements  of  the  size  of  lesions  in  B-mode  and  strain  images  and  the 
ratio  of  their  size  (strain/B-mode) 


Lesion  type 

Cyst 

Fibroadenoma 

IDC 

B-mode 

Width  (mm) 

9.8  ±  5.3 

11.8  ±4.1 

11.7  ±7.2 

Height  (mm) 

7.2  ±  3.4 

8.4  ±  2.4 

9.2  ±  4.8 

Area  (mm^) 

69.3  ±  72.2 

81.9  ±  44.4 

104  ±  135 

Strain 

Width  (mm) 

9.8  ±5.1 

11.6  ±4.1 

18.6  ±  7.3 

Height  (mm) 

7.6  ±  3.8 

8.8  ±  2.7 

14.8  ±  6.7 

rea  (mm^) 

72.8  ±  77.8 

81.2  ±  44.6 

236  ±  206 

Ratio  (strain/B-mode) 

Width 

1.02  ±0.16 

1.00  ±0.18 

1.74  ±  0.36 

Height 

1.04  ±0.16 

1.07  ±  0.21 

1.68  ±  0.36 

Area 

1.05  ±  0.24 

1.00  ±0.17 

3.02  ±  1.19 

Tabulated  values  are  the  mean  ±  the  SD  of  the  group.  The  size  ratio  for  benign  lesions  is  near  unity,  showing 
that  these  lesions  typically  have  the  same  size  in  both  imaging  modalities.  However,  invasive  ductal  carcinomas 
(IDC)  typically  are  2  to  3  times  as  large  in  strain  images  than  in  B-mode  images. 


was  necessary  to  observe  a  sequence  of  B-mode  and 
strain  images  displayed  side-by-side.  With  that  sequence, 
a  very  reproducible  determination  of  the  lesion  boundary 
could  be  obtained.  Measurements  of  lesion  dimension 
were  then  made  and  the  results  for  lesion  width  are 


0>) 

Fig.  10.  B-mode  and  strain  images  of  a  cluster  of  breast  cysts 
with  their  perimeter  traced  in  the  strain  image;  that  tracing  also 
in  the  B-mode  image  for  comparison. 


consistent  with  those  reported  by  Garra  et  al.  (1997). 
That  report  stated  a  lack  of  confidence  in  their  measure¬ 
ments  of  lesion  height.  Our  results  with  spherical  targets 
in  phantoms  show  that  we  can  accurately  measure  lesion 
dimension  in  both  height  and  width  (Zhu  and  Hall  2002) 
and,  therefore,  we  use  lesion  area  as  the  criterion  for 
comparing  lesion  size  in  B-mode  and  strain  images. 

The  significant,  but  monotonic,  change  in  strain 
contrast  as  a  function  of  precompression  appears  to  be 
unique  to  fibroadenomas  so  far  in  our  experience.  This 
contrast  variation  suggests  that  the  stress-strain  relation¬ 
ship  for  fibroadenoma  does  not  parallel  that  of  the  sur¬ 
rounding  tissue.  Fibroadenomas  that  vary  in  strain  con¬ 
trast  appear  dark  (stiffer)  at  low  precompression  and  lose 
contrast  (become  relatively  softer)  at  higher  precompres¬ 
sion.  One  possible  explanation  is  that  the  stress-strain 
relationship  for  the  surrounding  tissue  is  more  nonlinear 
than  that  of  the  fibroadenoma  over  the  range  that  each  are 
compressed  in  this  technique.  The  average  age  of  the 
women  with  fibroadenomas  in  this  study  was  44  years, 
and  their  dominant  breast  tissue  type  was  subjectively 
judged  to  be  fibroglandular  from  B-mode  images.  The 
data  reported  by  Krouskop  et  al.  (1998),  demonstrate  that 
both  glandular  tissue  and  (primarily)  fibrous  tissue,  such 
as  fibroadenoma,  have  nonlinear  stress-strain  relation¬ 
ships.  When  preloaded  to  5%  strain,  fibrous  tissue  is 
about  3  times  more  stiff  than  glandular  tissue.  However, 
the  appropriate  comparison  for  our  application  is  when 
both  tissues  have  minimal  preload  and  when  glandular 
tissue  is  preloaded  about  15%  and  fibrous  tissue  is  pre- 
loaded  some  small  fraction  of  that  (because  the  less  stiff 
tissues  strain  more  when  they  are  treated  as  a  composite), 
and  that  composite  is  strained  an  average  of  about  1%  for 
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data  acquisition.  Krouskop  and  colleagues  did  not  report 
those  results. 

Garra  et  al.  (1997)  tested  numerous  criteria  for  dis¬ 
criminating  benign  from  malignant  breast  lesions.  The 
criterion  that  provided  the  greatest  discrimination  in  their 
study  was  the  comparison  of  lesion  width  measured  in 
B-mode  and  strain  images.  They  attributed  that  size 
difference  to  the  desmoplasia  that  surrounds  most  ma¬ 
lignant  breast  tumors  (Tavassoli  1999).  Desmoplasia  is 
the  excessive  growth  of  fibrous  connective  tissue  in  the 
stroma  surrounding  the  malignancy.  That  growth  appears 
gray-white  and  feels  very  hard  in  gross  pathology 
(Tavassoli  1999).  Our  study  tested  this  criterion  and 
extended  the  observation  to  a  comparison  of  lesion  area. 
The  sequence  of  B-mode  and  strain  image  pairs  allows 
the  sonographer  to  select  images  representative  of  the 
“typical”  strain  image  for  a  lesion.  This  ability,  along 
with  better  determination  of  lesion  boundary  available  by 
viewing  a  sequence  of  images,  has  likely  improved  the 
ability  to  measure  true  lesion  size  in  strain  imaging 
compared  with  the  results  reported  by  Garra  et  al.  (1997). 
A  study  to  estimate  the  intraobserver  and  interobserver 
variability  in  choosing  the  “typical”  strain  image  and 
measuring  lesion  size  is  underway.  That  study  is  an 
essential  part  of  determining  the  value  of  the  relative  size 
of  lesions  in  B-mode  and  strain  image  pairs  as  a  diag¬ 
nostic  criterion.  The  utility  of  elasticity  imaging  in  dif¬ 
ferentiating  (from  benign  growths)  malignancies  that 
lack  desmoplasia  has  not  been  tested. 

Garra  et  al.  (1997)  also  found  the  brightness  of  the 
lesion  in  strain  images  to  be  a  useful  parameter,  but  our 
observation  of  the  changing  contrast  with  compression  in 
fibroadenomas  provides  an  improved  description  of  the 
contrast  of  solid  lesions  in  strain  imaging.  The  change  in 
strain  image  contrast  with  applied  compression  (e.g..  Fig. 
3)  demonstrates  that  observing  only  a  single  B-mode  and 
strain  image  pair  at  an  unknown  precompression  could 
be  very  misleading  in  the  interpretation  of  the  strain 
image  data.  A  sequence  of  side-by-side  B-mode  and 
strain  image  pairs  greatly  adds  to  the  ability  to  interpret 
strain  images. 

Numerous  other  criteria  were  tested  by  Garra  et  al. 
(1997),  but  demonstrated  limited  utility.  As  with  their 
study,  the  current  study  is  limited  by  the  small  number  of 
patients  and  lesions  included,  and  by  the  fact  only  one 
observer  was  involved  in  each  report.  A  larger  cohort  of 
patients  and  a  larger  number  of  observers  are  needed  to 


improve  the  statistical  analysis  of  fiiis  technique.  That 
effort  will  be  the  subject  of  a  future  report. 

CONCLUSIONS 

A  new  system  for  real-time  imaging  of  tissue  strain 
in  vivo  using  freehand  scanning  is  described  and  some  of 
the  results  obtained  with  this  system  are  reported.  The 
new  system  provides  real-time  feedback,  allowing  the 
user  to  manipulate  the  conditions  of  tissue  compression 
resulting  in  the  ability  to  successfully  scan  all  patients  for 
which  the  technique  was  attempted.  The  strain  images 
for  various  lesion  types  are  unique,  and  the  relative  size 
of  the  lesions  appears  to  be  a  useful  criterion  for  dis¬ 
criminating  benign  from  canceroiis  lesions.  However, 
further  testing  will  be  needed  to  support  this  observation. 
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This  manuscript  reports  a  technical  innovation  that  has  been  developed  for  real-time,  freehand  strain 
imaging.  This  work  is  based  on  a  well-known  block-matching  algorithm  with  two  significant  modifica¬ 
tions.  First,  since  displacements  are  estimated  row-by-row,  displacement  estimates  from  the  previous 
row  are  used  to  predict  the  displacement  estimates  in  the  current  row,  thereby  drastically  reducing  the 
search-region  size  and  increasing  computational  efficiency.  Second,  a  displacement  error  detection  and 
correction  method  is  developed  to  overcome  the  local  tracking  errors  that  may  be  more  severe  with  free¬ 
hand  scanning  and  thereby  improve  the  robustness  of  the  algorithm.  This  algorithm  has  been  imple¬ 
mented  on  a  clinical  ultrasound  imaging  system,  and  with  real-time  imaging  feedback,  long  sequences 
of  high  quality  strain  images  are  observed  using  freehand  compression.  Displacement  estimation  errors 
with  this  method  are  experimentally  measured  and  compared  with  results  from  simulation.  We  report 
only  a  specific  implementation,  with  no  comparison  to  other  displacement  estimators  in  the  literature 
and  no  optimization  of  this  specific  technique.  Images  of  tissue-mimicking  phantoms  with  small  spheri¬ 
cal  targets  are  used  to  test  the  ability  to  detect  small  lesions  using  the  strain  imaging  technique.  In  vivo 
strain  images  of  breast  and  thyroid  are  also  shown. 

Key  Words:  Elasticity;  elastography;  palpation;  tissue  characterization;  ultrasound. 


1.  INTRODUCTION 

Ultrasonic  strain  imaging^'’ is  expected  to  have  great  potential  for  improving  soft  tissue  di¬ 
agnosis.  The  vast  majority  of  the  strain  imaging  work  in  the  literature  has  focused  on  proof 
of  concept  and  algorithm  development.  However,  two  major  advances  need  to  occur  to 
make  ultrasonic  strain  imaging  a  clinically  useful  tool.  First,  there  is  a  need  to  acquire  rf  (or 
equivalent)  echo  data  under  clinically  acceptable  conditions  with  high  patient  throughput 
and  low  rescanning  rates.  Second,  there  is  a  need  to  further  develop  the  clinical  interpreta¬ 
tion  and  significance  of  these  results.  In  this  work,  we  concentrate  our  effort  on  the  first  is¬ 
sue. 

A  skilled  sonographer  can  acquire  B-mode  images  with  relative  ease.  Low  noise  strain 
images  are  more  difficult  to  produce.  Deformation  of  heterogeneous  tissue  can  result  in 
complex  motion  and  echo  signal  decorrelation  between  thepre-  and  postcompression  rf  echo 
frame  pair.  Echo  signal  decorrelation  leads  to  large  strain  estimation  errors.  Most  strain  im¬ 
aging  techniques  in  the  literature  produce  strain  images  in  two  steps.  First,  the  rf  echo  data  is 
acquired  by  either  digitizing  the  output  from  a  modified  clinical  ultrasound  system  or  di¬ 
rectly  downloading  the  rf  data  from  the  system  if  it  is  available.  Second,  the  stored  rf  echo 
frames  are  processed  off-line  to  produce  strain  images.  With  off-line  processing,  it  is  diffi¬ 
cult  to  know,  while  scanning,  whether  the  acquired  echo  frame  pairs  are  coherent  enough  for 
strain  estimation.  Hence,  off-line  data  processing  is  not  as  efficient  for  clinical  utilization  as 
real-time  imaging. 
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A  more  clinically  desirable  system  would  be  capable  of  displaying  B-mode  images  and 
strain  images  side-by-side  in  real-time.  Real-time  strain  imaging  can  help  clinicians  deter¬ 
mine  whether  satisfactory  strain  images  are  being  acquired  and  can  also  help  clinicians  im¬ 
prove  their  freehand  compression  techniques.  With  real-time  feedback,  clinicians  can  adjust 
the  compression  speed  and  angle  to  compensate  for  irregular  motion  (large-scale  rotation  or 
lateral  or  elevational  translation).  Hence,  real-time  strain  imaging  increases  data  acquisition 
efficiency. 

Another  important  aspect  for  clinical  utilization  of  strain  imaging  is  the  technique  used  to 
deform  the  soft  tissue.  Most  phantom  experiments  in  the  literature  used  motorized  compres¬ 
sion  fixtures.  These  devices  are  cumbersome,  limit  the  locations  that  strain  imaging  can  be 
applied,  and  are  time  consuming  to  incorporate.*  Thus  freehand  scanning  is  desirable.  Doy- 
ley  et  al  have  shown  that,  compared  with  motorized  compression,  the  penalty  for  freehand 
scanning  is  small  with  tissue  mimicking  phantoms.®  Hiltawsky  et  al  have  shown  that  strain 
imaging  of  breast  tissues  with  freehand  compression  is  feasible.’ 

Lorenz  et  al  have  developed  a  real-time  freehand  strain  imaging  system  that  is  external  to 
the  ultrasound  imaging  platform.’  Some  encouraging  results  are  obtained  for  prostate  appli¬ 
cations  using  the  system.  In  their  system,  a  modified  1-D  cross-correlation  algorithm  is  used 
for  displacement  estimation.  It  performs  well  with  small  strains  (much  less  than  1%)  since 
the  motion  is  tracked  only  in  one  dimension.  It  is  well-known  in  the  literature  that  the  con- 
trast-to-noise  ratio  in  strain  images  increases  with  the  applied  strain  (below  about  5%).’’*® 
Hence,  larger  single  step  compression  is  desirable. 

In  our  work,  strain-imaging  software  is  implemented  within  a  high-end  commercial  ultra¬ 
sonic  imaging  platform  (SONOLINE  Elegra,  Siemens  Medical  Solutions,  Ultrasound 
Group).  The  strain  images  are  displayed  in  real-time,  side-by-side  with  the  regular  B-mode 
images.  For  the  purpose  of  estimating  displacement  with  relatively  large  applied  strain 
{l'-2%  for  in  vivo  tissue  and  up  to  5%  for  tissue  mimicking  phantoms),  a  modified  2-D  block 
matching  algorithm  was  selected.  Block  matching  (template  matching)  algorithms  are 
widely  used  in  image  processing  applications  for  tracking  motion.  The  most  notable  appli¬ 
cation  is  video  compression  standards  such  as  MPEG.  Its  utilization  in  ultrasonic  imaging 
was  first  reported  by  Trahey  et  al  for  blood  flow  estimation.** 

The  purpose  of  this  paper  is  to  report  our  work  in  developing  a  real-time  freehand  strain 
imaging  technique.  We  report  only  a  specific  implementation  with  no  comparison  to  other 
displacement  estimators  in  the  literature  and  no  optimization  of  this  specific  technique. 

In  the  next  section,  we  will  introduce  our  displacement  and  strain  estimation  algorithms. 
In  the  results  section,  we  provide  performance  measures  of  this  algorithm  in  the  form  of  dis¬ 
placement  estimation  error  and  spatial  resolution  with  tissue-mimicking  phantoms.  We  also 
provide  examples  of  in  vivo  strain  images  obtained  from  this  system.  The  conclusion  section 
summarizes  this  work. 


METHODS 


1.  Standard  block  matching 

The  2-D  block  matching  algorithm  computes  a  sum-squared  difference  (SSD)  between 
pre-  and  pdstcompression  rf  frames  for  a  rectangular  kernel  over  a  search  region  as  follows: 

{K^/2  (K^n 

SSD(u,v)=  ^  ['i(/  +  i,J  +  y)-r2(/  +  /  +  M,y  +  y  +  v)]^ 


(1) 
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FIG.  1  Illustration  of  a  kernel  and  search  region,  (a)  Illustration  of  a  kernel,  (b)  The  search  region  defining  the 
range  of  locations  of  kernel  centers.  Each  ‘x’  represents  an  rf  data  sample.  A  kernel  is  composed  of  several  adjacent 
rf  A-lines  (5  in  our  implementation)  with  several  rf  samples  (1 1  in  our  implementation)  per  line  in  the  pre- 
compression  rf  field.  The  search  region  defines  the  range  of  possible  locations  of  the  center  of  the  kernel  in  the 
postcompression  rf  field. 


where  and  Tj  are  pre-  and  postcompression  rf  echo  fields,  respectively;  J  and  Jaic  axial  and 
lateral  if  sample  indices  for  the  location  where  the  displacement  is  estimated;  u  and  v  define 
search  locations  in  a  search  region;  and  and  are  kernel  height  and  width,  respectively. 
The  kernel  size  is  empirically  chosen  to  be  1 1  axial  samples  by  5  A-lines  (X'j  =  1 1  and  ATj  =  5) 
for  a  7.5MHz  transducer,  or  approximately  half  the  area  of  the  pulse-echo  point-spread  func¬ 
tion.  A  pictorial  illustration  of  the  kernel  and  search  region  is  shown  in  figure  1 .  For  each  lo¬ 
cation  (I,J)  at  which  the  displacement  is  estimated,  the  SSD  function  is  computed  for  every  rf 
sample  location  that  is  within  the  search  region  (range  of  kernel  centers)  defined  by 
-Uf=u=U2  and  -Fj=v=F2,  where  C/j,  ^2  represent  search  up,  down,  left  and  right 

distances,  respectively,  as  shown  in  figure  1.  The  search-region  height  and  width  are 
and  F=Fj+F2+l,  respectively.  The  displacement  distribution  usually  does  not 
need  to  be  estimated  as  finely  as  rf  samples.  We  use  A:  and  /  as  axial  and  lateral  indices  of  dis¬ 
placement  estimates.  The  location,  d2{Kl)\  at  which  the  minimum  SSD  is 

found  is  considered  to  be  the  displaced  position  of  the  kernel.  Hence,  and 
axial  and  lateral  displacement  estimates,  respectively. 

The  computational  cost  of  the  block-matching  algorithm  to  produce  one  displacement  es¬ 
timate  is  mainly  determined  by  the  kernel  and  the  search-region  sizes.  In  fact,  to  estimate 
one  displacement  vector,  the  subtraction-square-accumulation  operation  defined  in  Eq.  (1) 
needs  to  be  performed  times.  The  computed  SSD  values  are  dien  compared  UV 

times  to  find  the  minimum.  Since  the  kernel  size  is  usually  predefined  and  fixed,  the  task  of 
reducing  the  computational  cost  of the  block  matching  algorithm  is  to  find  away  to  minimize 
(/and  V, 

It  is  straightforward  to  estimate  the  computational  load  of  the  typical  implementation  of 
SSD-based  block  matching.  In  this  case,  U and  V are  selected  to  be  sufficiently  large  to  guar¬ 
antee  that  the  displacement  vector  is  enclosed  by  the  search  region.  Assume  the  size  of  each 
rf  echo  field  is  40  mm  by  40  mm  (typical  in  our  experiments)  and  the  maximum  applied  strain 
is  5%.  The  associated  axial  displacement  magnitude  is  then  0.05  x  40  mm = 2  mm.  The  axial 
sampling  frequency  for  our  study  is  36  MHz,  so  the  maximum  axial  displacement  magnitude 
is  about  94  samples.  Assuming  uniaxial  compression  of  an  incompressible  medium,  the  to- 
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tal  lateral  strain  is  at  most  5%  (assuming  no  elevational  expansion).  However,  it  is  difficult 
to  accurately  estimate  the  maximum  lateral  displacement  magnitude  since  we  do  not  know 
the  location  of  the  center  of  the  compression  (location  where  lateral  displacement  is  zero).  If 
we  assume  the  center  of  the  compression  occurs  at  the  middle  of  the  data  field,  the  maximum 
lateral  displacement  magnitude  is  0.05  x  20  mm  =  1  mm.  Given  that  the  rf  echo  field  usually 
consists  of  200  A-lines,  the  maximum  lateral  displacement  magnitude  is  about  5  A-lines. 
With  these  assumptions,  the  search-region  size  can  be  chosen  as  94x2+1  by  5x2+1  or  UV= 
189x11=2079. 

2.  Search  region  reduction 

The  size  of  the  search  region  can  be  minimized  by  using  prior  knowledge.  The  axial  and 
lateral  strain  are  defined  by  the  following  partial  differential  equations 


5i(^,/)  = 


^2(^.0  = 


dd,{k,l) 

dxi  ’ 

dd2(k,i) 


X2 


(2) 


In  the  sampled  space,  the  axial  and  lateral  strain  can  usually  be  approximated  by  the  fol¬ 
lowing  difference  equations 


(jt  n  d,{k,l)-d,(k-\,l)  (3) 


where  x^{kj)  and  X2(kJ)  are  the  axial  and  lateral  coordinates  of  the  location  of  the  displace¬ 
ment  estimate  respectively;  ^-land  M  represent  the  indices  of  adjacent  displacement 
estimates  relative  to  k  and  L  The  following  inequalities  can  be  derived  from  Eq.  (3) 


where  and  are  the  maximum  (allowed)  local  strain  magnitudes  in  the  axial  and  lateral  di¬ 
rections,  respectively.  The  spatial  separations  of  displacement  estimates  zxQXy{k,l)-x^{k-\J) 
and  X2(^,/)-X2(^,/-l),  respectively.  For  real-time  strain  imaging,  the  axial  separation  is  16 
samples  and  the  lateral  is  2  A-lines.  Experiments  in  phantoms  have  demonstrated  that  dis¬ 
placement  can  be  estimated  when  the  applied  strain  is  greater  than  5%.  However,  strain  con- 
trast-to-noise  ratio  decreases  rapidly  for  applied  strain  in  excess  of  5%,^^  Our  in  vivo 
experiments  have  shown  that  the  displacement  can  be  successfully  estimated  for  applied 
strain  up  to  about  2%.  Noise  dominates  in  the  displacement  estimates  when  the  applied 
strain  is  more  than  2%  for  in  vivo  tissues.  We  can  use  this  as  prior  knowledge  to  limit  die  ex- 
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ROI 

FIG.  2  Illustration  of  an  ROI.  Grid  points  (circles)  arc  locations  at  which  the  displacement  distribution  is  esti¬ 
mated.  Each  grid  point  coincides  with  an  rf  sample  location,  but  the  displacement  distribution  is  less  densely  sam¬ 
pled  than  the  rf  echo  signals.  In  our  real-time  implementation,  grid  point  separations  are  16  samples  in  the  axial  di¬ 
rection  and  2  A-lines  in  the  lateral  direction. 


tent  of  the  search  region.  Note  that  if  the  displacement  difference  between  adjacent  esti¬ 
mates  is  1  sample  in  the  axial  direction  and  1  A-line  in  the  lateral  direction,  then  Inequalities 
(4)  give  us 


Si  <  6.25%,  (5) 

5'2<50%. 


The  maximum  local  strain  we  intend  to  estimate  is  less  than  62SVo.  If  we  estimate  dis¬ 
placement  in  the  order  of  increasing  k  and  /  (row  by  row  and  from  left  to  right),  then  the  dis¬ 
placement  that  is  currently  being  estimated  is  within  a  3  sample  by  3  sample  block  centered 
at  the  location  predicted  by  the  previously  estimated  displacement  (one  sample  in  each  direc¬ 
tion  away  from  our  best  guess).  In  other  words,  we  can  reduce  the  search-region  size  to  3  by 
3  ifwe  use  the  previous  estimates  to  predict  where  to  search  for  adjacent  estimates.  Note  that 
we  allow  larger  local  strain  than  the  total  strain  since  the  strain  distribution  is  not  uniform. 
With  the  search  region  reduction,  I7F-3x3=9.  Compared  to  typical  block  matching,  the  new 
method  reduces  the  computational  load  by  a  factor  of 2079/9  =  231. 

In  implementing  the  reduced  search-region  block-matching  strategy,  we  first  manually  se¬ 
lect  a  region  of  interest  (ROI)  which  is  a  subregion  of  the  field  of  view.  For  example,  in 
breast  imaging,  we  select  an  ROI  that  excludes  undesired  echo  regions.  The  locations  at 
which  the  displacement  is  estimated  are  determined  by  grid  points  with  equal  spacing  start¬ 
ing  at  the  upper-left  comer  of  the  ROI,  as  shown  in  figure  2.  The  displacement  is  then  esti¬ 
mated  in  two  stages.  In  the  first  stage,  the  displacement  of  the  first  row  of  grid  points,  as 
shown  in  figure  2,  is  estimated  using  Eq.  (1).  Since  there  is  no  prior  knowledge  of  the  dis¬ 
placement  distribution,  a  large  search  region  is  used.  The  size  of  the  search  region  at  this 
stage  is  determined  by  the  following  equations. 


Ui^U2=^CEIL{SiXi{Qfi)\ 

Vi=V2=^CElL{SiA\ 


(6) 
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where  CEIL  is  the  function  that  rounds  to  the  next  larger  integer;  Xj(0,0)  is  the  depth  of  the  top 
of  the  ROI;  A  is  the  number  of  A-lines  in  the  field  of  view.  The  search  region  created  by  Eq. 
(6)  is  large  enough  to  enclose  the  true  displacements  as  long  as  the  applied  strain  does  not  ex¬ 
ceed  iS'j. 

In  the  second  stage,  the  displacement  is  estimated  from  the  top  of  the  ROI  to  the  bottom, 
row  by  row.  Displacement  estimates  obtained  in  the  first  row  are  used  to  predict  displace¬ 
ments  in  lower  rows,  and  the  search  region  in  this  stage  is  reduced  to  3  if  samples  by  3  A-lines 
using  the  equation 

SSD(u,v)=  Y,  Ylr,{I+i,J+j)- 

r^il  +  d^ik-l,l)  +  i  +  u,J +d2ik-\,l)+ j +v)f . 


In  Eq.  (7),  a  search  center,  (/+</, (A:-l,  /)?  /))>  is  used  to  guide  the  search.  In  other 

words,  when  we  do  not  have  any  knowledge  of  the  displacement  distribution,  the  search  cen¬ 
ter  is  (/,  J),  After  a  row  (or  rows)  of  the  displacement  distribution  is  estimated,  that  informa¬ 
tion  can  be  used  to  guide  the  search  at  neighboring  locations  and  allows  the  use  of  a  small 
search  region.  Note  from  Inequalities  4  that  it  is  logical  to  set  the  search-region  center  to  be 
{Rd^{k-\^  /),  J^djik,  /-!)).  However,  since  the  difference  between  dj(^kA,  1)  and  is 

small,  the  search-region  center  is  selected  as  {I+d^{kA,  /),  J+d^{kA,  t))  to  simplify  the  algo¬ 
rithm. 

The  computational  load  can  be  further  reduced  by  performing  the  2-D  (3  by  3  search  re¬ 
gion)  search  sparsely.  Lateral  displacement  does  not  need  to  be  estimated  as  densely  as  axial 
displacement  because  lateral  sample  spacing  is  much  greater  than  axial.  For  real-time  imag¬ 
ing,  we  can  apply  3  by  3  search  regions  every  5  rows  and  use  the  lateral  displacement  esti¬ 
mates  to  guide  the  next  4  rows  of  displacement  estimation.  These  4  rows  of  displacement 
estimates  are  obtained  using  a  3  by  1  search  region  (1-D  search).  An  even  more  aggressive 
strategy  is  to  only  estiinate  the  lateral  displacement  for  the  first  row  and  use  this  to  predict  the 
lateral  displacement  for  the  remaining  rows  and  limit  the  2-D  search  to  the  first  row  only. 

3.  Displacement  error  detection  and  correction 

Two  types  of  displacement  estimation  errors  can  occur.  It  is  common  for  large  errors  to 
appear  in  the  first  row.  This  is  due  to  the  inherent  pre-  and  postcompression  rf  waveform 
decorrelation  and  periodic  ambiguities  (‘false-peak  errors’)  associated  with  large  search  re¬ 
gions.  Figure  3  shows  an  example  of  this  type  of  error. 

The  second  type  of  error  results  from  correlations  in  displacement  estimates  when  using 
predicted  displacements  to  reduce  the  search-region  size.  Errors  in  displacement  estimates 
propagate  if  they  are  large  enough  that  the  defined  search  regions  do  not  enclose  true  dis¬ 
placements.  When  tissue  is  compressed,  large  and  irregular  local  deformation  can  occur. 
This  may  cause  local  decorrelation  in  the  recorded  rf  frame  pair.  Figure  4  shows  an  example 
of  this  type  of  error.  The  displacement  estimation  errors  tend  to  accumulate  when  estimating 
displacement  near  those  local  regions. 

A  segmentation  method  that  can  detect  and  correct  large  errors  is  needed  to  overcome 
these  problems.  Each  row  of  the  displacement  estimates  is  checked  for  errors  in  three  steps. 
In  the  first  step,  from  left  to  right  within  a  row,  the  displacement  estimates  are  segmented. 
Segmentation  occurs  if  the  difference  between  adjacent  displacement  estimates  is  larger 
than  2  samples.  The  result  of  this  step  for  the  displacement  curve  shown  in  figure  3  (row 
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FIG.  3  Example  of  errors  in  the  first  row  displacement  estimates. 


Lateral  Distance  (mm) 

FIG,  4  Example  of  errors  in  consecutive  rows  with  a  3x3  search  region. 


/I+4)  is  shown  in  figure  5(a)  where  the  displacement  estimates  are  segmented  into  three 
groups.  In  the  second  step,  groups  that  are  not  adjacent  are  merged  if  the  difference  of  dis¬ 
placement  estimates  between  two  nearest  end  points  of  the  two  groups  is  smaller  than  a 
threshold  (3  samples  in  this  case).  For  the  example  displacement  curve,  group  3  is  merged 
into  group  1 .  In  the  last  step,  the  group  that  has  the  largest  number  of  displacement  estimates 
(‘members’,  group  1  in  this  example)  is  marked  as  the  ‘correct’  group  of  displacement  estimates 
and  all  remaining  groups  are  marked  as  errors.  In  the  error  correction  stage,  the  displacement 
values  of  the  error  groups  are  then  discarded  and  replaced  by  linearly  interpolating  values 
fi-om  the  correct  group.  Figure  5(b)  shows  the  displacement  curve  after  error  correction. 
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FIG.  5  Demonstration  of  error  detection  and  correction,  (a)  Groups  that  are  generated  after  initial  segmentation. 
This  displacement  curve  is  the  same  as  row  n+4  in  figure  4.  (b)  Displacement  curve  after  error  correction  with  the 
segmentation  method. 


Tests  with  in  vivo  data  have  shown  that  the  error  detection  and  correction  process  does  not 
need  to  be  applied  to  eveiy  row  of  displacement  estimates.  Displacement  errors  build  up 
gradually  since  the  small  search  region  prevents  large  displacement  deviations  between  ad¬ 
jacent  rows.  Errors  are  more  apparent  and  easier  to  detect  after  the  displacement  estimation 
process  progresses  several  rows.  Widi  this  observation,  we  apply  the  error  detection  and  cor¬ 
rection  process  once  every  5  rows.  For  each  detected  error,  all  5  displacement  estimates  are 
replaced  by  the  interpolated  values. 

4.  Subsample  accuracy  displacement  estimation 

The  displacement  distribution  that  is  estimated  using  this  modified  block  matching  algo¬ 
rithm  has  integer  sample  accuracy.  With  36  MHz  sampling  and  the  strain  estimation  method 
described  below,  we  find  in  phantom  experiments  that  when  the  total  applied  strain  is  larger 
than  2%,  this  accuracy  is  adequate  for  creating  low  noise  strain  images.  However,  when  the 
total  applied  strain  is  smaller  than  2%,  obvious  strain  artifacts  can  be  seen  in  the  image.*^” 
There  are  two  ways  of  alleviating  diis  problem.  One  way  is  to  interpolate  recorded  if  frames 
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PIG.  7  Photograph  of  a  transducer  with  the  compressor  plate  mounted. 


to  higher  sampling  frequencies.  This  will  increase  the  computational  cost  tremendously 
since  interpolation  requires  computation  and  the  kernel  and  search-region  sizes  must  be  in¬ 
creased.  The  alternative  method,  used  in  our  implementation,  is  to  quadratically  interpolate 
the  SSD  values  around  the  minimum  to  obtain  sub-sample  accuracy  in  displacement  esti¬ 
mates.  With  this  added  processing,  the  strain  artifacts  are  not  severe  when  the  applied  strain 
is  more  than  0.2%. 

5.  Strain  estimation 

The  axial  strain  is  defined  as  the  spatial  derivative  of  the  axial  displacement,  and  there  are 
several  methods  to  estimate  this  derivative.  The  simplest  methods  are  forward,  backward, 
and  center  differences  where  only  two  data  points  are  used.  Estimating  derivatives  using 
only  two  data  points  requires  low  noise  in  displacement  estimates.  Since  a  relatively  small 
kernel  is  used  to  estimate  displacement,  the  noise  in  the  displacement  is  too  high  to  use  these 
methods.  However,  axial  strain  can  be  estimated  using  a  low  order  polynomial  curve  fitting 
method,^^  and  we  have  implemented  a  linear  regression  strain  estimator.  In  addition,  this 
method  provides  the  ability  to  trade  off  spatial  resolution  for  increased  smoothness  of  strain 
images.  The  strain  images  have  better  spatial  resolution,  but  more  noise,  if  shorter  segments 
of  displacement  estimates  are  used  in  linear  regression.  The  strain  images  are  visually  well 
balanced  in  smoothness  and  spatial  resolution  if  the  linear  window  length  is  set  between  2-3 
mm  for  a  7,5  MHz  transducer. 


IMPLEMENTATION 

We  have  implemented  the  modified  block  matching  algorithm  on  the  Siemens  Sonoline 
Elegra.  The  strain  imaging  software  is  an  application  that  resides  in  the  Elegra.  The 
real-time  beamformed  I-Q  (anal)1ical  representation  of  the  rf  signal)  frames  are  passed  to  a 
digital  signal  processor  subsystem.  That  subsystem  houses  two  Texas  Instruments  TMS- 
320C80  MVP  processors  tiiat  execute  the  software.  An  I-Q  frame  pair  is  used  for  displace¬ 
ment  and  strain  estimation.  The  first  frame  is  also  envelope  detected  and  aB-mode  image  is 
formed.  Then,  both  B-mode  and  strain  images  are  displayed  side-by-side  on  the  Elegra’s 
monitor  as  shown  in  figure  6.  The  strain  image  corresponds  to  a  region  of  interest  (ROI) 
marked  by  the  white-outlined  rectangular  subregion  on  the  B-mode  image.  A  user  interfacie 
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FIG.  8  Plots  of  the  variance  in  displacement  estimate  errors  for  experimental  and  simulated  data,  (a)  Motorized 
compression,  (b)  Freehand  compression,  (c)  Curve-fits  for  experimental  data  with  simulated  data.  In  (c),  the  solid 
and  dotted  curves  represent  results  obtained  fi'om  motorized  and  freehand  compression,  respectively;  circles  are 
data  obtained  from  simulation. 


allows  the  adjustment  of  the  size  and  location  of  the  ROI,  the  separation  (in  the  data  stream) 
between  I-Q  frames  used  to  estimate  strain,  and  some  of  the  strain  visualization  parameters 
such  as  strain  to  gray-scale  mapping,  etc.  This  software  is  capable  of  displaying  the 
side-by-side  images  at  about  7  frames/second. 

With  this  system,  strain  imaging  is  performed  in  three  stages.  In  the  setup  stage,  the 
sonographer  locates  the  lesion  and  selects  the  ROI  in  which  the  strain  is  estimated  and  dis¬ 
played.  Then,  by  pushing  a  button,  the  software  enters  the  real-time  side-by-side  display 
mode.  The  user  starts  the  compress-release  cycle.  If  the  user  finds  an  image  sequence  that  is 
of  interest,  she  can  freeze  the  data  acquisition.  A  cine-mode  allows  the  user  to  browse  each 
frame  or  loop  through  a  selected  set  of  frames  for  a  more  careful  study  of  the  acquired  data. 

Tissue  deformation  is  generated  by  cyclic  motion  of  the  transducer  (i.e.,  compressing  and 
releasing  the  body  surface).  A  small  (4.5  cm  x  9  cm)  compressor  plate  can  be  mounted  to  the 
face  of  the  transducer  to  enlarge  the  compression  surface  and  produce  a  more  uniform  stress 
field^^  A  picture  of  the  transducer  with  a  mounted  compressor  plate  is  shown  in  figure  7. 
Krouskop  et  al^’  have  shown  that  if  the  cyclic  motion  is  approximately  IHz,  then  the  breast 
tissue  components  behave  as  elastic  materials  (i.e.,  the  viscous  effects  are  negligible). 
Hence,  during  data  acquisition,  cyclic  deformation  of  about  IHz  is  attempted. 


REAL-TIME  FREEHAND  STRAIN  IMAGING 


11 


Experiments  have  shown  that  7/s  frame  rate  is  sufficient  for  a  sonographer  to  scan  free¬ 
hand  and  control  the  compression  motion  of  the  transducer  to  compensate  for  imdesirable 
lateral  and  elevational  motion.  No  additional  motion  restricting  device  is  necessary.  How¬ 
ever,  successfully  acquiring  strain  image  data  is  not  trivial.  The  sonographer  needs  to  under¬ 
stand  that  the  information  being  extracted  is  mechanical  in  nature  and  that  the  major 
challenge  in  strain  image  scanning  is  to  minimize  the  rf  waveform  decorrelation. 

Note  that  most  of  the  parameters  involved  in  the  processing  that  generates  strain  estimates 
have  not  been  optimized.  These  parameters  include  the  kernel  size  used  in  the  block¬ 
matching  algorithm  and  the  window  length  used  in  the  moving  linear  regression  algorithm 
that  estimates  strain  from  the  displacement  distribution.  The  optimization  will  be  performed 
in  our  future  work. 


RESULTS 

All  data  sets  shown  in  this  section  were  acquired  using  the  system  described  above.  A 
7.5L40  linear  array  transducer  (pulsed  at  7.2MHz)  was  used  in  our  data  acquisition.  The 
field  of  view  was  40  mm  x  40  mm.  The  system  performance  is  studied  in  terms  of  basic  im¬ 
age  quality  parameters.  Reported  here  are  representative  measurements  of  noise  and  resolu¬ 
tion.  More  detailed  investigations  into  these  topics  will  follow  in  future  work  after  the 
processing  parameters  are  optimized. 

1«  Displacement  estimation  error 

The  variance  of  the  error  in  the  displacement  estimates  was  measured  using  a  uniform  gel¬ 
atin  phantom  to  produce  a  predictable  displacement  distribution.  Both  motorized  compres¬ 
sion  and  freehand  compression  were  used  in  order  to  compare  the  variance  of  displacement 
errors  with  each  of  these  methods.  The  motorized  compression  (a  laboratory  system*®)  used 
a  large  compression  plate  that  covered  the  entire  upper  surface  of  the  free-standing  gelatin 
block(10cmx  1 0  cm  x  7  cm  (WxDxH),  no  additional  fixtures).  The  motor  was  programmed 
to  produce  a  sinusoidal  compression  of  20%  at  0.4  Hz.  Freehand  scanning  was  performed 
with  the  compressor  plate  shown  in  figure  7  and  the  compression  was  intended  to  replicate 
that  of  motorized  compression.  The  separation  between  I-Q  frame  pairs,  called  the  skip 
number  and  used  for  strain  image  formation,  was  adjusted  to  achieve  a  wide  range  of 
frame-average  strains  from  these  data  sets. 

The  acquired  data  sets  were  processed  off-line  using  an  algorithm  identical  to  that  pro¬ 
grammed  on  the  Elegra.  Since  the  phantom  had  uniform  stiffness,  the  displacement  curve 
along  the  compression  direction  should  be  a  straight  line  (strain  is  constant  over  the  entire 
field  of  view).  Linear  regression  was  applied  to  the  estimated  displacement  curve  along  each 
A-line  in  the  region  of  interest  to  generate  the  best-fit  displacement  curve.  This  line  was  then 
considered  to  be  the  true  displacement.  The  displacement  error  was  calculated  as  the  differ¬ 
ence  between  estimates  and  the  fitted  lines.  The  corresponding  strain  was  calculated  by  av¬ 
eraging  over  strain  estimates  for  all  A-lines  in  the  ROL 

We  also  simulated  rf  frame  pairs  for  the  medium  with  imiform  stiffness.^  The  scanning 
pulse,  sampled  at  36  MHz,  had  7,2  MHz  center  frequency,  a  -6  dB  bandwidth  of  40%  and  a 
Gaussian  lateral  profile  with  -6  dB  beam-width  of  400  pm  and  beam  spacing  of  200  pm. 
These  parameters  closely  simulated  the  beam  profile  produced  by  the  Elegra  7.5L40  probe. 
The  simulated  compressions  resulted  in  applied  strains  of  0.2, 0.5, 1, 2, 3, 4  and  5%.  For  each 
compression,  30  rf  frame  pairs  were  generated  and  the  modified  block-matching  algorithm 


(a)  (b) 

FIG.  9  Strain  image  (a)  of  apbantom  with  spherical  targets.  The  average  strain  is  2.5%.  The  size  measurements, 
shown  in  (b),  are  accurate,  as  detailed  in  table  1. 

was  used  for  displacement  estimation.  The  displacement  estimates  were  compared  with  the 
known  true  displacements  to  calculate  the  variance  of  the  displacement  estimation  error. 

Figure  8  shows  the  variance  of  the  displacement  estimation  error  versus  the  estimated  ap¬ 
plied  strain.  In  figure  8(a)  and  (b)  (motorized  and  fi*eehand  compression,  respectively)  the 
first  50  frames  in  the  collected  rf  data  sequences  were  used  as  precompression  data  fields. 
The  skip  number  was  varied  from  1  to  15  for  motorized  compression  and  from  1  to  25  for 
fi*eehand  compression  to  achieve  strains  ranging  from  0. 1%  to  5%.  There  were  a  total  of  725 
strain  and  variance  measurements  for  motorized  compression  and  1250  measurements  for 
freehand  compression.  In  figure  8(a),  there  are  6  measurements  with  high  displacement  esti¬ 
mation  error  at  relatively  high  strain.  These  are  die  cases  where  the  error  detection  and  cor¬ 
rection  method  failed  due  to  excessive  noise  in  the  first  row  displacement  estimates. 

Displacement  error  variance  estimates  for  motorized  and  freehand  compression  were  fit  to 
a  second  degree  polynomial  in  log-log  space  to  generate  representative  curves  for  each  data 
set  and  those  curves  were  plotted  in  figure  8(c).  Note  that  the  6  measurements  for  motorized 
compression  with  high  error  variance  were  excluded  when  curve  fitting.  The  circles  in  fig¬ 
ure  8(c)  are  displacement  error  variance  measurements  obtained  from  the  simulation.  The 
standard  deviations  of  the  error  variance  measurements  for  simulation  are  so  small  that  they 
are  not  plotted  (they  would  not  be  visible  if  plotted). 

As  seen  in  figure  8(c),  the  displacement  error  variance  curve  is  relatively  flat  for  strain  less 
than  1%.  Thisislikelyduetothefixeddisplacementerrorproducedbythequadratic interpo¬ 
lation.*^  The  displacenjent  estimate  error  variance  for  experimental  data  increases  with  the 
applied  strain  faster  flian  the  results  obtained  from  simulation.  This  is  likely  due  to 
elevational  motion  resulting  from  compressing  the  free-standing  gel  block  (plane  stress  con¬ 
ditions),  whereas  the  simulated  data  employed  plane  strain  conditions  (no  elevational  mo¬ 
tion).  Comparable  performance  is  observed  with  motorized  and  freehand  compression. 
Although  motorized  compression  generally  has  slightly  lower  displacement  errors,  the  ben¬ 
efit  associated  with  freehand  scanning  offsets  the  small  improvement  in  displacement  esti¬ 
mate  errors. 

2.  SmaU  lesion  detection 

Gelatin  phantoms  with  spherical  targets  that  were  three  times  stiffer  than  the  background** 
were  used  to  test  the  strain  imaging  system  performance  with  small  targets.  The  strain  image 
(acquired  with  2.5%  compression)  shown  in  figure  9(a)  contains  three  targets  (4.0  mm,  3.2 
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TABLE  1  Measured  target  size  in  millimeters. 


Leftmost  target 

Middle  target 

Rightmost  target 

Height  (axial) 

4.1 

3.1 

2.3 

Width  (lateral) 

4.0 

3.4 

2.3 

Actual  diameter 

4.0 

3.2 

2.4 

mm  and  2.4  mm  diameter,  respectively).  The  apparent  size  of  these  targets  was  measured  in 
the  strain  image  in  both  axial  and  lateral  directions.  Figure  9(b)  shows  line  segments  that 
correspond  to  the  width  and  height  for  each  target  and  table  1  shows  the  measured  sizes. 
These  results  suggest  that  spherical  targets  as  small  as  2.4  mm  diameter  can  be  accurately 
measured  in  both  the  lateral  and  axial  dimensions. 

3.  In  vivo  strain  images 

Real-time  freehand  strain  imaging  has  also  been  performed  on  in  vivo  tissues.  The  images, 
shown  below,  demonstrate  that  these  strain  images  have  reasonably  low  noise  and  high  con¬ 
trast. 

Figure  10  shows  an  in  vivo  breast  cyst  that  is  about  3  mm  x  3  mm.  The  visibility  of  the  cyst 
in  the  strain  image  confirms  our  phantom  results  that  lesions  of  a  few  millimeters  in  diameter 
can  be  detected  in  the  strain  images.  The  exact  reason  that  the  fluid  filled  cyst  appears  stiffer 
than  the  background  is  unknown.  A  reasonable  hypothesis  is  that  the  cyst  fluid  is  bounded  by 
a  distended  capsule  and  appears  stiff  much  like  an  air-filled  balloon  feels  stiff. 

Figure  1 1  shows  an  in  vivo  breast  carcinoma.  The  apparent  size  of  the  tumor  is  much  larger 
in  the  strain  image  (about  twice  as  big)  1h^  in  the  B-mode  image.  This  is  consistent  with  the 
findings  of  Garra  et  al.® 

Figure  12  shows  an  in  vivo  thyroid  strain  image.  There  is  a  nodule  inside  the  thyroid,  seen 
in  both  B-mode  and  strain  images.  The  tissue  structures  and  therefore  the  boundary  condi¬ 
tions  around  the  thyroid  are  very  different  from  the  breast,  and  in  both  cases  compression  in¬ 
duced  motion  is  complex.  However,  with  real-time  feedback,  the  sonographer  can  manipulate 
the  compression  technique  and  obtain  strain  images.  A  problem  in  strain  imaging  of  the  thy¬ 
roid  is  that  the  trachea  and  major  blood  vessels  are  often  included  in  the  field  of  view.  Since 


14 


ZHU  AND  HALL 


k 

»i  -.v-  4  ► 

^  L-  * 

■  . . :.U:. 

;  1  .  .. 

'^art^cj  .  Nodule  Tractiea 

Carotid  Nodule  Trachea 

FIG.  12  Strain  image  of  an  in  vivo  thyroid  with  a  nodule.  Flow  in  the  carotid  and  echo  noise  in  die  trachea  cause 
errors  in  motion  tracking  and  strain  estimation. 


there  are  no  echo  signals  from  the  trachea,  displacement  estimates  in  diis  region  are  at  best 
misleading.  The  blood  flow  in  the  carotid  is  perpendicular  to  the  image  plane  and  introduces 
elevational  motion  that  causes  echo  signal  decorrelation  and  motion  tracking  errors.  The  ob¬ 
server  must  consider  these  factors  for  correct  strain  image  interpretation  in  this  case. 


CONCLUSION 

A  computationally  efficient  displacement  estimation  algorithm  has  been  developed  for 
real-time,  freehand  ultrasonic  strain  imaging.  The  proposed  method  is  based  on  a  block¬ 
matching  algorithm  that  is  widely  used  for  motion  detection  in  digital  image  processing.  Ma¬ 
jor  modifications  increase  the  computational  efficiency  and  robustness  of  the  typical  block 
matching  algorithm.  The  algorithm  is  implemented  on  the  Siemens  SONOLINE  Elegra  as 
an  add-on  software  application. 

With  real-time  feedback  of  strain  images,  sonographers  can  adjust  their  compression/ 
scanning  technique  to  consistently  form  strain  images.  Strain  images  with  acceptable  qual¬ 
ity  are  observed  in  both  breast  and  thyroid  scanning,  which  require  different  scanning  tech- 
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niques  and  where  the  motion  is  much  more  complex  than  in  phantoms.  Since  the  algorithm  is 
implemented  on  a  commercially  available  clinical  imaging  system,  data  can  be  efficiently 
acquired  from  a  large  number  of  patients,  enabling  clinical  evaluation  of  strain  imaging  in 
soft  tissue  diagnosis. 


ACKNOWLEDGEMENTS 

We  are  grateful  for  the  technical  assistant  provided  by  Siemens  Medical  Solutions  Ultra^ 
sound  Group.  Without  their  help  this  work  would  not  have  been  accomplished.  All  in  vivo 
scanning  was  done  by  Candace  S,  Spalding.  We  greatly  appreciate  her  help.  We  are  also 
grateful  for  the  financial  support  of  USAMRAA  DAMD 17-00- 1-0596  and  NSF  BES- 
970822 1 .  The  U.S.  Army  Medical  Research  Acquisition  Activity,  820  Chandler  Street,  Fort 
Detrick  MD  21 702-5014  is  the  awarding  and  administering  acquisition  office  forDAMD17- 
00- 1 -0596 .  The  information  reported  here  does  not  necessarily  reflect  the  position  or  policy 
of  the  U.S.  Government,  and  no  official  endorsement  should  be  inferred. 


REFERENCES 

1 .  Tristam,  M.,  Barbosa,  D.C.,  Cosgrove,  D.O.,Nassiri,  D.K.,  Bamber,  LC.  and  Hill,  C.R,,  Ultrasonic  study  of  in 
vivo  kinetic  characteristics  of  human  tissues.  Ultrasound  Med.  Biol  12, 927-937  (1986). 

2.  Krouskop,  T.A.,  Dougherty,  D.R.  and  Vinson,  S.F.,  A  pulsed  Doppler  ultrasonic  system  for  making 
non-invasive  measurements  of  the  mechanical  properties  of  soft  tissues,  J.  Rehab.  Res.  Dev.  24, 1-8  (1987). 

3.  Ophir,  J.,  C6spedes,  I.,  Ponnekanti,  H.,  Vazdi,  Y.,  and  Li,  X.,  Elastography:  a  quantitative  method  for  imaging 
the  elasticity  of  biological  tissues.  Ultrasonic  Imaging  13,  1 1 1-134  (1991). 

4.  O’Donnell,  M.,  Skovoroda,  A.R,,  Shapo,  B,,  and  Emelianov,  S. ,  Internal  displacement  and  strain  iamging  us¬ 
ing  ultrasonic  speckle  tracking,  /EEE  Trans.  Ultrason.  Ferroelec.  Freq.  Contr.  "41, 314-325  (1994). 

5.  Chaturvedi,  P.,  Insana,  M.F.  and  Hall,  T.  J. ,  2-d  companding  for  noise  reduction  in  strain  imaging,  IEEE  Trans. 
Ultrason.  Ferroelec.  Freq.  Contr.  45,  \79-\9l  (199%). 

6.  Doyley,  M.M.,  Bamber,  J.C.,  Fuechsel,  F.  and  Bush,  N.L,,  A  freehand  elastogr^hic  imaging  approach  for 
clinical  breast  imaging:  system  development  and  performance  evaluation.  Ultrasound  Med.  Biol  27, 1347-1357 
(2001). 

7.  Hiltawsky,  K..,  Kruger,  M.,  Starke,  C,,  Heuser,  L.,  Ermert,  H.,  and  Jensen,  A,  Freehand  ultrasound 
elastography  of  breast  lesions:  clinical  results.  Ultrasound  Med.  Biol  27, 1461-1469  (2001), 

8.  Garra,  B.S.,  Cespedes,  I.,  Ophir,  J.,  Spratt,  S.R.,  Zuuibier,  R.A.,  Magnant,  C.M.  and  Pennanen,  M.F., 
Elastography  of  breast  lesions:  initial  clinical  results.  Radiology  202, 79-86  (1997). 

9.  Lorenz,  A.,  Sommerfeld,  H.  J.,  Garcia-Schurmann,  M.,  Philippou,  S.,  Senge,  T.  and  Ermert,  H. ,  A  new  system 
for  the  acquistioin  of  ultrasonic  multicompression  strain  images  of  the  human  prostate  in  vivo,  IEEE  Trans. 
Ultrason.  Ferroelec.  Freq.  Contr.  46,1147-1154(1999). 

1 0.  Chaturvedi,  P.,  Insana,  M.F,,  and  Hall,  T.J.,  Testing  the  limitations  of  2-D  local  companding  in  strain  imaging 
using  phantoms, /EjEE  Trans.  Ultrason.  Ferroelec.  Freq.  Contr.  45,  1022^1 03 1  (1998). 

11.  Trahey,  GE.,  Allison,  J.  W.  and  von  Ramm,  O.T.,  Angle  independent  ultrasonic  detection  of  blood  flow, 
IEEE  Trans.  Biomed  Eng.  34, 965-967  (1987). 

12.  Varghese,  T.  and  Ophir,  J.,  Characterization  of  elastogr^hic  noise  using  the  envelope  of  echo  signals.  Ultra¬ 
sound  Med  Biol  24,543-555  (1998). 

1 3 .  Alam,  S.K.  and  Ophir,  J. ,  The  effects  of  nonlinear  signal  transforamtions  on  bias  errors  in  elastography,  IEEE 
Trans.  Ultrason.  Ferroelec.  Freq.  Contr.  47, 297-303  (2000). 

14.  Kallel,  F.  and  Ophir,  J.,  A  least-squares  strain  estimator  for  elastography,  Ultrasonic  Imaging  10,  195-208 

(1997).  , 

15.  Urkowitz,  H.,  Signal  Theory  and  Random  Processes  (Artech  House,  Inc.,  Norwood,  MA,  1983). 


16 


ZHU  AND  HALL 


16.  Konofagou,  E.,  Dutta,  P.,  Ophir,  J.  and  C6spedes,  I.,  Reduction  of  stress  nonuniformities  by  apodization  of 
compressor  displacement  in  elastography.  Ultrasound  Med.  Biol  22^  1229-1236  (1996). 

17.  Krouskop,  T.A.,  Wheeler,  T.M.,  Kallel,  F.,  Garra,  B.S.  and  Hall,  T.J.,  Elastic  moduli  of  breast  and  prostate  tis¬ 
sues  under  compression.  Ultrasonic  Imaging  20, 260-274  (1998). 

18.  Hall,  T.J.,  Bilgen,  M.,  Insana,  M.F.,  and  Krouskop,  T.A.,  Phantom  materials  for  elastography,  IEEE  Trans 
Ultrason,  Ferroelec.  Freq  Contr.  44, 1355-1365  (1997). 


890 


IEEE  TRANSACTIONS  ON  MEDICAL  IMAGING,  VOL.  22,  NO.  7,  JULY  2003 


A  Finite-Element  Approach  for  Young’s  Modulus 
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Abstract — ^Modulus  imaging  has  great  potential  in  soft-tissue 
characterization  since  it  reveals  intrinsic  mechanical  properties. 
A  novel  Young’s  modulus  reconstruction  algorithm  that  is  based 
on  finite-element  analysis  is  reported  here.  This  new  method  over¬ 
comes  some  limitations  in  other  Young’s  modulus  reconstruction 
methods.  Specifically,  it  relaxes  the  force  boundary  condition  re¬ 
quirements  so  that  only  the  force  distribution  at  the  compression 
surface  is  necessary,  thus  making  the  new  method  more  practical. 
The  validity  of  the  new  method  is  demonstrated  and  the  perfor¬ 
mance  of  the  algorithm  with  noise  in  the  input  data  is  tested  using 
numerical  simulations.  Details  of  how  to  apply  this  method  under 
clinical  conditions  is  also  discussed. 

Index  Terms — Tissue  characterization,  tissue  elasticity. 


I.  Introduction 

The  elastic  properties  of  biological  tissues  are  usually 
modified  by  disease.  Surgeons  often  describe  the  “  feel” 
of  excised  abnormal  tissues.  As  a  result,  a  quantitative  mea¬ 
sure  of  the  elastic  properties  of  tissue  should  be  useful  in  di¬ 
agnosing  abnormalities.  The  physical  quantities  that  describe 
tissue  elastic  properties  are  stress,  strain,  and  elastic  moduli,  and 
methods  have  been  developed  to  estimate  each  of  these.  Palpa¬ 
tion,  which  has  been  used  for  more  than  4000  years,  utilizes 
tissue  surface  stress  information  to  detect  tissue  abnormalities. 
Palpation  remains  an  effective  diagnostic  tool.  In  fact,  the  ma¬ 
jority  of  breast  tumors  are  discovered  with  palpation  [1].  How¬ 
ever,  palpation  is  qualitative  and  lacks  sensitivity  to  small  deep 
abnormalities.  Quantitative  methods  similar  to  palpation  have 
been  developed  to  visualize  surface  pressure  [2],  [3].  Other  re¬ 
cent  developments  in  bioelasticity  imaging  techniques  involve 
accurately  and  noninvasively  measuring  the  tissue  strain  distri¬ 
bution  during  external  compression.  Studies  have  shown  that 
these  techniques  show  promise  in  diagnosing  and  monitoring 
diseases  of  the  breast  i4]-[7],  kidney  [8]-[  11],  and  blood  vessels 
[12],  [13]. 

Mapping  stress  or  strain  distributions  provides  only  relative 
information  about  tissue  elasticity.  Using  either  stress  or  strain 
information  alone,  one  can  only  identify  a  region  of  tissue  that  is 
stiff  (or  soft)  relative  to  its  surroundings.  Elastic  moduli  provide 
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an  absolute  measure  of  tissue  elasticity  that  is  intrinsic  to  the  ma¬ 
terial.  The  stress  or  strain  distributions  alone  lacks  a  one-to-one 
relationship  with  the  elastic  moduli  distribution.  Images  of  the 
stress  or  strain  distribution  may  also  include  misleading  arti¬ 
facts  that  could  lead  to  uncertainties  in  diagnosing  tissue  abnor¬ 
malities.  Therefore,  it  is  desirable  to  measure  elastic  moduli  in 
bioelasticity  imaging  techniques.  However,  measuring  the  dis¬ 
tribution  of  elastic  moduli  is  more  difficult  than  either  the  stress 
or  strain  distribution. 

The  theory  of  mechanics  shows  that  to  describe  the  com¬ 
plete  elastic  properties  of  a  material  requires  a  tensor  that  has 
81  components  [14],  Clearly,  it  is  impractical  to  measure  all 
these  components.  Assumptions  can  be  made  to  simplify  the 
problem  and  reduce  the  number  of  unique  tensor  elements.  If 
a  material  is  assumed  to  be  continuous,  incompressible,  and 
isotropic,  then  its  elasticity  can  be  completely  described  by  one 
elastic  modulus,  either  Young’s  modulus  E  or  shear  modulus  p,. 
Strictly  speaking,  none  of  the  above  assumptions  are  valid  for 
biological  tissues,  but  most  biological  tissues  closely  approx¬ 
imate  continuous  and  incompressible  materials.  Some  tissues, 
such  as  muscle,  are  anisotropic  in  their  structure,  function,  and 
mechanical  properties.  For  this  paper,  however,  we  will  assume 
tissue  to  be  continuous,  incompressible,  and  isotropic  as  a  first 
approximation. 

Currently,  ultrasonic-based  techniques  for  measuring  the 
elastic  modulus  of  tissue  fall  into  two  categories.  First,  dynamic 
compression  techniques  [I^HIS],  such  as  sonoelasticity,  use 
a  vibrator  to  propagate  low-ffequency  “  pumping”  waves  into 
tissue.  In  the  most  promising  of  these  approaches,  shear  wave 
velocity  or  wavelength  are  estimated,  and  from  these  the  shear 
modulus  can  be  estimated.  However,  problems  associated  with 
this  technique  are  high  image  noise,  low  spatial  resolution,  and 
difficulty  in  propagating  the  shear  wave  eneigy  across  tissue 
boundaries. 

The  other  category  is  referred  to  as  (quasi)static  compres¬ 
sion  techniques.  In  static  compression  techniques,  the  tissue 
Young’s  modulus  distribution  is  estimated  from  the  tissue  de¬ 
formation  and  boundary  pressure  measurements.  The  methods 
to  estimate  tissue  deformation  have  been  extensively  discussed 
in  ultrasound  based  elastography  [19]-[28],  The  tissue  is 
deformed  either  by  an  external  force  or  an  internal  force.  The 
RF  echo  waveforms  before  and  after  an  incremental  deforma¬ 
tion  are  recorded,  and  the  tissue  displacement  distribution  is 
estimated  by  comparing  these  RF  waveforms.  Tissue  internal 
displacement  can  be  also  obtained  using  magnetic  resonance 
imaging  [29]-[31]  and  optical  elastography  [32]  techniques. 
Young’s  modulus  estimation  can  be  performed  utilizing  the 
tissue  deformation  information  obtained  with  the  strain  imaging 
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techniques.  In  addition  to  the  displacement  distribution,  some 
Yoxmg’s  modulus  estimation  methods  also  require  knowledge 
of  the  pressure  or  force  boundary  conditions. 

There  are  four  methods  in  the  literature  for  reconstructing  the 
Yoxmg’s  modulus  distribution  based  on  static  compression  tech¬ 
niques  for  displacement  estimation.  The  first  method  estimates 
Young’s  modulus  by  numerically  solving  a  second-order  partial 
differential  equation  that  describes  a  linear,  isotropic,  incom¬ 
pressible  medium  under  static  deformation  [33].  That  method 
requires  significant  spatial  smoothing  of  the  displacement  esti¬ 
mates  to  obtain  second-order  partial  differentials  that  are  also 
smooth.  Hence,  with  noisy  displacement  estimates,  that  method 
inherently  has  low  spatial  resolution.  Another  problem  asso¬ 
ciated  with  that  method  is  that  for  a  two-dimensional  (2-D) 
analysis,  the  force  boundary  condition  of  the  medium  must  be 
known  on  all  sides.  However,  in  practice,  the  force  distribution 
can  only  be  (easily)  measured  on  one  side  (the  compression  sur¬ 
face)  of  the  medium. 

The  second  method  uses  an  iterative  technique  to  reconstruct 
the  modulus  distribution  [34],  [35].  That  method  uses  finite- 
element  analysis  (FEA)  to  solve  the  forward  elasticity  problem. 
The  input  to  the  FEA  algorithm  is  the  measured  displacement 
field,  the  assumed  boundary  conditions,  and  an  initial  guess  of 
the  modulus  distribution.  The  output  of  the  FEA  algorithm  is 
an  estimate  of  the  displacement  distribution.  The  difference  be¬ 
tween  the  measured  displacement  distribution  and  the  FEA  pre¬ 
diction  is  used  to  adjust  the  modulus  distribution  from  its  initial 
guess.  By  repeating  the  process  multiple  times,  one  can  obtain 
a  modulus  distribution  that  minimizes  the  displacement  distri¬ 
bution  difference  in  a  least  squares  sense.  The  advantage  of  that 
approach  is  that  it  does  not  require  knowledge  of  the  pressure 
boundary  conditions.  However,  without  knowing  the  boundary 
pressure,  only  relative  modulus  estimates  can  be  obtained.  In 
other  words,  the  ratio  of  the  modulus  between  different  loca¬ 
tions  can  be  determined.  Although  that  method  can  reduce  the 
artifacts  in  strain  images,  it  does  not  provide  absolute  measure¬ 
ment  of  the  tissue  modulus  distribution  which  can  be  useful  in 
tumor  discrimination  as  suggested  in  [36],  and  an  incorrect  ini¬ 
tial  modulus  guess  may  result  in  convergence  to  an  incorrect 
modulus  distribution.  For  media,  such  as  tissue,  that  have  a  com¬ 
plicated  modulus  distribution,  a  good  initial  guess  for  the  mod¬ 
ulus  distribution  is  difficult  to  obtain. 

In  the  third  modulus  reconstruction  method,  a  finite-differ¬ 
ence  approach  is  used  to  describe  the  elasticity  problem  in  a 
medium  [37].  That  approach  rearranges  linear  equations  that 
describe  the  forward  problem  so  that  the  modulus  distribution 
becomes  unknown  variables  in  these  equations.  The  modulus 
distribution  can  then  be  solved.  However,  that  method  also  re¬ 
quires  knowledge  of  the  boundary  conditions  on  all  sides  of  the 
object. 

The  fourth  approach  uses  a  variational  method  to  formulate 
the  forward  solution  [38].  Then  the  terms  with  unknowns  are 
rearranged  to  derive  a  matrix  equation  similar  to  ours.  However, 
the  boundary  force  condition  was  not  utilized  in  their  treatment. 
Hence,  this  method  can  only  reconstruct  the  ratio  between  the 
Lame  constants  and  tissue  mass  density. 

In  our  approach,  FEA  is  used  to  construct  a  set  of  linear  equa¬ 
tions  that  describes  the  elastic  behavior  of  a  2-D  object.  Similar 


to  the  third  method  mentioned  above,  we  rewrite  the  linear  equa¬ 
tion  set  so  that  the  Young’s  modulus  distribution  are  explicit 
variables  which  can  be  solved.  The  solution  does  not  require 
an  initial  guess  or  iteration  of  the  modulus  distribution  solution, 
and  it  provides  absolute,  not  relative,  modulus  estimates.  Unlike 
the  equation  set  for  solving  forward  elasticity  problems,  where 
the  number  of  equations  equals  the  number  of  unknown  vari¬ 
ables,  the  equation  set  for  our  inverse  solution  usually  involves 
more  equations  than  unknowns.  This  allows  us  to  simplify  the 
force  boundary  conditions  so  that  only  one  (surface)  force  dis¬ 
tribution  is  necessaiy  to  solve  for  the  modulus  distribution. 

The  details  of  our  modulus  estimation  method  are  described 
in  Section  II.  The  validity  of  this  method  is  tested  with  simu¬ 
lations  and  results  are  shown  in  Section  HI.  The  discussion  of 
how  this  technique  can  be  implemented  for  ultrasonic  imaging 
systems  is  provided  in  Section  IV. 

II.  Method 

Three  integral  parts  of  the  proposed  modulus  estimation 
method  are  described  in  this  section,  one  subsection  each. 
Section  II-A  provides  the  information  necessary  for  solving 
a  forward  elasticity  problem  of  2-D  continua  using  FEA. 
Although  the  content  of  this  subsection  is  well  known  in 
the  literature,  it  is  briefly  reviewed  here  to  provide  sufficient 
background,  terminology,  and  notation  for  the  development 
of  Section  II-B.  In  Section  II-B,  the  FEA-based  modulus 
estimation  technique  (inverse  problem)  is  described  in  detail. 
Section  II-C  addresses  issues  of  how  to  apply  the  proposed 
method  under  practical  constraints. 

A,  The  FEA  Method  for  Solving  a  Forward  Problem  of  2-D 
Continua 

The  FEA  procedure  for  solving  a  forward  elasticity  problem 
of  2-D  continua  can  be  summarized  as  follows  [39]-[41]. 

1)  Select  an  element  type  and  derive  the  element  stiffness 
matrix. 

2)  Form  a  mesh  using  the  selected  element  to  cover  the  re¬ 
gion  of  interest  (ROI)  for  which  the  elasticity  problem  is 
solved. 

3)  Generate  the  global  stiffness  matrix  by  assembling  ele¬ 
ment  stiffness  matrices. 

4)  Apply  the  boundary  conditions  to  solve  the  global  matrix 
equations  for  the  solution. 

The  details  of  these  steps  are  provided  below. 

Step  1:  For  2-D  problems,  the  common  choices  for  element 
type  are  triangles  or  quadrilaterals.  In  other  words,  each  element 
has  either  three  or  four  nodes.  The  element  matrix  equations  for 
elasticity  problems  have  the  form 

/^(e)5(e)  ^  fie) 

where  and  are,  respectively,  the  element  stiffness 

matrix,  the  element  nodal  displacement  vector,  and  the  element 
nodal  force  vector  for  element  e.  In  this  paper,  rectangular 
elements,  as  shown  in  Fig.  1 ,  are  used.  Details  of  how  to  compute 
the  element  stiffhess  matrix  is  provided  in  Appendix  I. 

Step  2:  For  problems  that  can  be  described  by  partial  differ¬ 
ential  equations  but  do  not  have  closed  form  solutions,  FEA  has 
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Fig.  1.  The  rectangular  mesh  element.  Nodal  numbers  start  from  the  bottom 
left  comer  and  increase  on  the  clockwise  direction. 


been  developed  to  find  approximate  numerical  solutions  on  dis¬ 
cretized  problem  domains.  These  solutions  can  then  be  interpo¬ 
lated  to  form  continuous  solution  spaces  using  shape  functions. 
The  discretization  process  is  performed  by  creating  a  mesh  that 
covers  the  ROI  in  the  object. 

A  mesh  is  composed  of  elements  that  cover  a  contiguous  area 
in  the  problem  domain.  In  this  paper,  all  elements  in  the  ROI  are 
rectangles  of  the  same  size.  A  nine-element  mesh  is  illustrated 
in  Fig.  2.  The  numbers  in  the  center  of  the  rectangular  elements 
are  element  numbers  and  numbers  close  to  nodes  (intersection 
points)  are  nodal  numbers. 

The  mesh  configuration  can  be  represented  by  a  connectivity 
matrix,  (7.  The  number  of  rows  of  C  equals  the  number  of  el¬ 
ements,  Ne-  For  four-node  elements  (rectangles),  C  has  four 
columns.  The  ith  row  of  C  records  nodes  associated  with  the 
ith  element.  For  the  mesh  shown  in  Fig.  2,  the  first  two  rows  of 
C  are 

ar,  O' 

Step  5;  The  matrix  equation  for  a  meshed  system  (for  ex¬ 
ample,  the  system  described  in  Fig.  2)  has  the  form 

K8  =  /,  (3) 

where  K  is  the  global  stifihess  matrix;  8  is  the  global  nodal 
displacement  vector;  and  /  is  the  global  nodal  force  vector.  The 
global  displacement  vector  has  the  form 

^  ( <^lx  8iy  82x  ^2y  •  •  -  (4) 

where  81^  and  8iy  are  the  displacements  in  the  x  and  y  directions 
for  node  1 ,  and  so  on.  The  global  nodal  force  vector  has  the  form 

f  =  {fix  fly  hx  f2y  **  *  (5) 

where  fix  and  fiy  are  the  net  force  exerted  on  node  1  in  x  and  y 
directions,  and  so  on.  The  global  stiffness  matrix  K  is  assembled 
fi*om  element  stiffness  matrices.  The  assembly  process  can  be 
found  in  Appendix  II. 

Step  4:  Displacement  boundary  conditions  are  usually  speci¬ 
fied  for  the  elasticity  problems  in  our  applications,  and  a  penalty 
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Fig.  2.  A  nine-element  mesh  composed  of  rectangular  elements  with  uniform 
size.  Numbers  in  the  center  of  the  rectangular  elements  are  the  element  numbers 
and  numbers  close  to  nodes  (intersection  points)  are  the  global  nodal  numbers. 

approach  [4 1  ]  can  be  used  to  solve  (3).  The  details  of  the  penalty 
approach  is  provided  in  Appendix  III. 

B,  The  FEA  Approach  for  Solving  an  Inverse  Problem  of  2-D 
Continua 

Soft  tissues  can  generally  be  considered  as  incompressible 
media  [33].  Hence,  the  Possion’s  ratio  can  be  assumed  to  be  a 
constant  that  is  close  to  0.5  (0.49,  for  instance)  throughout  the 
ROI.  With  this  assumption,  matrix  K  in  (24)  is  same  for  every 
element  if  all  elements  have  the  same  aspect  ratio  (which  is  true 
in  this  paper  since  the  mesh  is  composed  of  the  elements  with 
same  size). 

In  the  element-to-global  stiffness  matrix  assembly  procedure 
described  in  Section  11- A,  each  component  of  (i.e.,  the 
product  of  the  element  Young’s  modulus  and  a  constant)  is  ac¬ 
cumulated  onto  the  global  stiffness  matrix.  Hence,  each  compo¬ 
nent  of  the  global  stiffness  matrix  is  a  linear  combination  of  the 
Young’s  modulus  of  each  element.  In  other  words,  the  global 
stiffness  matrix  can  be  written  as 

K  =  [Kij\  *  =  =  (6) 

iv. 

(7) 

e-1 

where  are  constants. 

From  (3),  the  left-hand  side  of  the  system  matrix  equation  is 
K8.  In  Young’s  modulus  reconstruction,  the  displacement  dis¬ 
tribution  is  estimated  with  tissue  motion  tracking  techniques. 
In  other  words,  ^  is  a  “known”  vector.  Performing  the  matrix- 
vector  multiplication 

•  (8) 
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Since  Kij  is  a  linear  combination  of  [see  (7)],  each  compo¬ 

nent  of  the  resulting  vector  on  the  right-hand  side  of  (8)  is  also 
a  linear  combination  of  .  Or 


K6^ 


(9) 


The  right-hand  side  of  (9)  can  be  written  as  a  product  of  a 
matrix  and  a  vector.  Or 


K8  =  DE  (10) 

where  D  is  a  iV-by-iVe  matrix; 
E  =  •••  •••  f  is  the  Young’s 

modulus  vector.  Now,  (3)  can  be  rewritten  as 

DE^f.  (11) 

Recall  from  the  global  stiffness  matrix  assembly  procedure, 
provided  in  Section  II-B,  that  the  component  K\f  = 
of  the  local  stiffness  matrix  of  element  e  is  accumulated  to  the 
Cl  '^th  row  and  ^th  column  of  the  global  stiffness  matrix,  or 

Kcie)  ^(e) .  When  performing  the  multiplication  ofKd  [see  (8)], 

E^^^Kij  is  multiplied  with  S^(e) .  The  product  is  then  accumu- 

lated  to  the  ^th  component  of  the  resulting  vector.  Hence, 
Kij8^(e)  is  a  summand  of  the  C^^^th  row  and  eth  column  com¬ 
ponent  of  the  matrix  D.  Based  on  this  observation,  the  matrix 
D  can  be  assembled  from  the  element  stiffness  matrix  using  the 
following  procedure. 

1)  Initialize  a  iV-by-iVe  null  matrix. 

2)  For  element  e,  generate  a  local  variable  number  to  global 
variable  number  conversion  index  vector  defined  in 
(26). 

3)  Accumulate  Kij8ji^)  to  for  i  =  1, . . . ,  8  and  j  = 

4)  Iterate  2  and  3  for  all  elements. 

Similar  to  solving  forward  elasticity  problems  with  (3),  the 
Young’s  modulus  reconstruction  problem  can  be  solved  from 
(1 1)  were  D  is  an  N-by~Ne  matrix.  Usually,  N  >  Ne  and  (1 1) 
defines  an  overdetermined  set  of  equations.  The  common  tech¬ 
nique  for  solving  an  overdetermined  linear  equation  set  is  to 
convert  it  to  a  least-square  problem  [42].  The  conversion  can 
be  done  by  multiplying  both  sides  of  (1 1)  by  D'^ 

{D‘^D)E  =  D'^f.  (12) 

Since  (D^D)  is  an  N^-hy-Ne  matrix,  E  can  be  solved  by  the 
following  equation 

E  =  (13) 

The  application  of  this  method  under  practical  constraints  is 
provided  in  the  next  subsection. 

C.  Practical  Concerns 

1)  Necessary  Measurements:  Equation  (3)  implies  that 
there  is  an  unstressed  state  to  which  the  object  returns  when  all 
external  forces  are  removed  from  the  object.  Then,  the  external 


forces  (e.g.,  gravitational  force,  atmospherical  pressure,  and 
compressional  force)  are  exerted  on  the  object.  The  displace¬ 
ment  is  measured  between  the  unstressed  state  and  the  state 
with  external  load.  In  reality,  however,  the  geometrical  distri¬ 
bution  of  the  object  in  the  unstressed  state  is  unknown.  Hence, 
(11)  [which  is  derived  from  (3)]  cannot  be  used  directly  to 
solve  the  inverse  problem.  Fortunately,  if  the  object  is  assumed 
to  be  a  linear  elastic  body  for  small  incremental  deformations, 
then  this  problem  can  be  solved. 

Assume  there  are  two  loading  states  of  the  object  Si  and  S2. 
In  iSi,  the  object  can  be  described  by 

(14) 

where  6^^'  ^  is  the  object  displacement  between  the  natural  state 
and5i;  is  the  compressional  load  measured  in  Si.  J11S2, 
the  object  can  be  described  by 

KS(S2)  ^  f(S2)^  (15) 

Note  that  K  is  state  independent  given  the  linear  elastic  body 
assumption. 

Subtracting  (15)  from  (14),  we  find 

KA8  3=  A/  (16) 

where  A(^  =  8^^'^^  —  and  A/  =  Replacing 

(3)  with  (16),  the  derivation  introduced  in  Section  II-B  still 
holds.  (16)  states  that  three  measurements,  namely,  A8, 
and  ,  are  necessary  for  Young’s  modulus  reconstruction. 

2)  Young's  Modulus  Reconstruction  With  Partial  Boundary 
Conditions:  Recall  that  (11)  defines  an  over-determined  set 
of  linear  equations.  This  means  that  only  a  subset  of  (1 1)  are 
needed  for  inverse  solution.  Let  Ds  be  a  matrix  that  is  formed 
from  a  subset  of  rows  from  D,  and  fs  be  the  force  vector 
formed  from  the  corresponding  subset  of  /.  As  long  as  (Dj’ Ds) 
is  invertible,  a  unique  inverse  solution  can  be  obtained  from 

E^iDfDsr^D-^U  (17) 

This  property  can  be  used  to  relax  the  force  boundary  conditions 
for  the  inverse  problem. 

Under  typical  conditions  a  subregion  of  the  tissue  imder  study 
is  observed.  A  mesh  can  be  created  to  cover  an  ROI  which 
is  a  subset  of  the  field  of  view.  Only  one  side  (the  surface  of 
the  tissue)  of  the  boundary  force  distribution  can  be  measured 
with  ease  with  a  force  sensor  array.  With  the  example  shown  in 
Fig.  2,  let  us  assume  that  the  compressional  force  can  be  mea¬ 
sured  along  the  top  side  (nodes  {13,  14,  15,  16}).  Ds  and  fs 
are  determined  by  the  following  rule:  select  all  equations  from 
(11)  involving  all  interior  nodes  of  the  mesh  and  all  nodes  ex¬ 
cept  the  outer  most  two  for  which  the  force  measurements  are 
made.  With  the  example  shown  in  Fig.  2,  the  selected  nodes  are 
(6,7, 10, 11, 14, 15}.  The  related  subset  ofthe  matrix  equations 
include  rows  {11,  12, 13, 14, 19,  20, 21,  22,  27,  28,  29,  30}  of 
(11). 

It  is  difficult  to  mathematically  prove  that  the  equation  se¬ 
lection  rule  described  above  always  provides  invertible  Ds, 
however,  it  is  easy  to  test  whether  D^  Ds  is  invertible  given  a 
mesh  configuration.  From  the  variety  of  mesh  configurations 
that  have  been  tested,  the  equation  selection  rule  described 
above  always  provides  invertible  D^Dg. 
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3)  The  Size  of  D  Matrix:  The  size  of  the  matrix  D  can  be 
a  problem  if  it  is  assembled  as  a  dense  matrix.  For  example, 
if  the  mesh  has  100  x  100  elements,  then  there  are  total  of 
101  X  101  =  10201  nodes.  The  size  of  D  is  20402-by-10000, 
or  204020  000  components.  Representing  each  component 
as  double  precision  floating  point  numbers  and  storing  this 
dense  matrix  requires  about  1.52  GB  of  computer  memory. 
With  current  computer  technology,  storing  such  a  large  matrix 
may  not  be  a  problem.  However,  calculating  is 

impractical. 

Fortunately,  H  is  a  sparse  matrix.  Recall  that  each  row  of  (3) 
describes  the  behavior  of  a  node,  and  there  are  a  maximum  of 
four  elements  related  to  one  node.  Thus,  the  maximum  nonzero 
components  of  a  row  in  D  is  four. 

The  number  of  nonzero  components  of  D  for  a  lOO-by-100 
element  mesh  is  in  fact  80  000.  Using  the  sparse  matrix  features 
provided  by  MATLAB,  this  can  be  stored  with  about  1MB  of 
memory.  Compared  with  dense  matrix  storage,  this  is  a  1600:1 
reduction  in  memory  requirement.  Using  a  750-MHz  Pentium 
m  PC,  can  be  computed  in  about  20  s.  Note  that 

although  the  discussion  is  based  on  D,  same  conclusion  can  be 
drawn  for  Ds. 


III.  Results 

A.  Solving  the  Forward  Problem  With  Ideal  Input 

First,  we  simulated  an  object  for  which  the  forward  problem 
was  solved.  The  dimension  of  the  simulated  object  was 
40  X  40  X  2  (width  x  height  x  thickness  in  millimeters). 
Young’s  modulus  distribution  of  the  object  is  shown  in  Fig.  3. 
The  Young’s  modulus  of  the  background  was  15  kPa  which 
approximates  the  stiffness  of  normal  glandular  breast  tissue 
[36].  There  were  two  10-mm-diameter  targets  in  the  object  that 
simulate  lesions.  The  upper  target  was  three  times  stiffer  (45 
kPa)  than  the  background  and  the  lower  target  was  three  times 
less  stiff  (5  kPa)  than  the  background. 

A  mesh  was  created  for  this  object  with  160  elements  in  both 
horizontal  and  vertical  directions.  The  total  number  of  elements 
was  160  X 160  =  25  600.  The  size  of  each  element  was  0.25  mm 
x  0.25  mm.  To  produce  a  realistic  estimate  of  boundary  force 
values,  we  assumed  plain  stress  conditions  for  the  compression. 
The  Possion’s  ratio  u  of  0.49  was  used  (incompressible  media). 
The  displacement  boundary  conditions  were  assigned  such  that 
the  displacement  of  bottom  side  of  the  object  was  zero  and  the 
top  side  of  the  the  object  was  0.8  mm  simulating  a  2%  compres¬ 
sion  of  the  object. 

Using  the  forward  FEA  method  introduced  in  Section  II-A, 
we  calculated  the  displacement  distribution  and  the  force  distri¬ 
bution  along  top  and  the  bottom  sides  of  the  object.  Using  the 
displacement  distribution,  strain  in  both  the  horizontal  direction 
(-Sxx)  and  vertical  direction  {syy)  were  calculated  and  shown  in 
Fig.  4(a)  and  (b),  respectively.  With  the  plain  stress  assumption, 
die  following  relationship  holds  s^xlsyy  —  -v  12,  The  force 
distribution  on  the  top  and  bottom  sides  of  the  object  are  shown 
in  Fig.  5(a)  and  (b),  respectively.  Note  that  the  forces  on  the  top 
side  are  all  negative  since  the  direction  of  the  applied  force  is 
pointing  vertically  down.  As  shown  in  Fig.  5,  the  magnitude  of 
force  on  the  outer  most  nodes  are  half  of  the  magnitude  of  the 


Fig,  3.  The  Young’s  modulus  distribution  for  the  simulated  object.  The  units 
of  the  color  bar  are  kPa, 


(b) 

Fig.  4.  Strain  images  obtained  from  the  forward  FEA  calculation  for  the  object 
illustrated  in  Fig.  3.  (a)  Horizontal  strain  (b)  Vertical  strain  Syy. 


force  on  their  adjacent  nodes  since  the  area  of  support  for  the 
outer  most  nodes  is  half  that  of  the  inner  nodes.  With  the  same 
surface  pressure,  the  exerted  force  is  half  the  magnitude. 
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(a) 


(b) 

Fig.  5.  The  boundary  force  distribution  obtained  from  the  forward  FEA 
calculation  for  tihe  object  illustrated  in  Fig.  3.  (a)  Force  on  the  top  side, 
(b)  Force  on  the  bottom  side. 


B.  Solving  the  Inverse  Problem  With  Ideal  Input 

To  test  the  modulus  reconstruction  technique  we  used  the  re¬ 
sults  of  the  forward  problem  calculations,  but,  in  effect,  dis¬ 
carded  nonessential  information.  The  only  input  to  the  inverse 
(modulus  distribution)  computation  were  the  ideal  displacement 
distribution  over  a  (sub-)ROI  and  the  force  distribution  at  the 
top  (compression)  surface  generated  by  the  forward  solution. 
The  sub-ROI  was  30  mm  x  30  mm  and  a  new  mesh  was  created 
to  cover  that  area.  An  illustration  of  the  meshed  areas  for  the 
forward  simulation  and  the  inverse  reconstruction  is  shown  in 
Fig.  6. 

Fig.  7  shows  the  result  of  the  Young’s  modulus  estimation. 
The  standard  deviation  of  the  relative  error  in  the  reconstructed 
Young’s  modulus  distribution  is  2.0  x  10”^% — ^a  nearly  exact 
reconstruction  is  obtained.  This  result  is  encouraging  since  it 
shows  that  the  Young’s  modulus  estimation  method  introduced 
above  is  a  valid  approach  and  confirms  that  the  “partial  force” 


Fig.  6.  An  illustration  of  the  meshed  areas  for  the  forward  simulation  and  the 
Young’s  modulus  reconstruction. 


Fig.  7.  The  reconstructed  Young’s  modulus  image  using  the  ideal  (noise-free) 
displacement  and  force  boundary  conditions. 

boundary  condition  is  sufficient  to  estimate  the  modulus 
distribution. 

In  the  above  modulus  estimation  simulation,  the  size  of  the 
elements  was  the  same  as  that  used  in  the  forward  simulation. 
The  small  elements  provide  high  spatial  resolution.  However,  it 
is  not  likely  that  such  high  spatial  resolution  can  be  achieved 
under  all  practical  conditions.  In  the  next  simulation,  we  in¬ 
creased  the  size  of  the  elements  to  1  mm  x  1  mm.  In  other 
words,  we  blurred  the  spatial  sampling  by  a  factor  of  4.  Fig.  8(a) 
shows  the  “true”  modulus  image  of  the  same  object.  The  “true” 
modulus  value  of  each  element  were  calculated  by  averaging  the 
16  modulus  values  in  the  corresponding  area  of  the  finer  meshed 
object.  The  inverse  problem  was  solved  using  the  exact  displace¬ 
ment  and  force  from  the  forward  simulation  results  produced  by 
the  finer  mesh.  The  result  of  the  modulus  estimation  is  shown 
in  Fig.  8(b).  The  relative  difference  between  Fig.  8(a)  and  (b) 
is  shown  in  Fig.  8(c).  The  mean  and  standard  deviation  of  the 
image  shown  in  Fig.  8(c)  are  —1.25%  and  17,8%,  respectively. 
The  small  mean  value  suggests  that  the  modulus  estimates  are 
unbiased.  However,  the  standard  deviation  value  shows  that  the 
reconstructed  modulus  image  using  larger  element  size  can  be 
noisy.  Since  the  modulus  estimates  are  unbiased,  a  simple  spa¬ 
tial  averaging  can  improve  the  visual  effect  of  the  reconstructed 
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(c) 

Fig.  8.  Simulation  results  using  the  larger  element  size,  (a)  The  ideal  decimated  modulus  image,  (b)  The  reconstructed  modulus  image  using  the  larger  element 
size,  (c)  The  relative  difference  between  (a)  and  (b)  measured  in  percentage,  mean  =  —1.25%,  std  =  17.8%. 


modulus  images  at  the  expense  of  further  reducing  the  spatial 
resolution. 

C  Modulus  Estimation  With  Noisy  Displacement  and 
Boundary  Force  Estimates 

Inevitably,  there  is  noise  in  the  displacement  and  boundary 
force  estimates  that  are  used  in  modulus  estimation.  Hence, 
it  is  necessary  to  study  the  effect  of  noise  in  the  input  data 
on  the  resulting  modulus  estimates.  To  avoid  the  difficulty 
of  analytically  deriving  the  noise  relationship  between  input 
data  and  final  outcome,  we  rely  on  numerical  simulation. 
The  forward  solution  shown  in  Figs.  4  and  5  was  used  as 
the  ideal  displacement  and  boundary  force  distributions.  The 
inverse  simulation  was  based  on  the  mesh  configuration  that 
produced  the  result  shown  in  Fig.  8.  Noise,  modeled  as  zero 
mean  white  Gaussian  random  processes  [43],^  was  added  to 

*Biigen  et  al  have  shown  through  simulation  that  the  noise  in  displacement 
estimates  is  Gaussian  distributed  However,  the  spectrum  of  the  noise  is 
not  shown  in  dieir  woiic.  The  spectrum  of  noise  is  assumed  to  be  white 
as  an  ^roximation. 


both  the  ideal  displacement  and  boundary  force  distributions,  A 
range  of  the  standard  deviations  (noise)  were  used  to  study  the 
relationship  between  noise  power  and  the  modulus  estimation 
error.  For  each  predetermined  level  of  noise,  100  realizations 
of  (noisy)  displacement  and  force  distributions  were  generated 
and  the  object  modulus  distributions  were  reconstructed. 

The  modulus  estimation  error  is  defined  as  the  difference  be¬ 
tween  the  “true”  modulus  distribution  [Fig.  8(a)]  and  the  esti¬ 
mated  modulus  distribution  Avith  noise  present  in  the  input  data 
(force  boundary  condition  and  displacement  distribution).  The 
relative  mean  and  the  relative  standard  deviation  of  the  error 
were  calculated  from  the  outcome  of  all  100  realizations  of  noise 
fields.  The  simulation  results  are  shown  in  Fig,  9.  The  standard 
deviation  of  the  noise  added  to  the  boundary  force  distribution 
is  5%  of  the  mean  ideal  force  in  Fig.  9(a)  and  (b)  and  10%  in 
Fig.  9(c)  and  (d).  Fig.  9(a)  and  (c)  shows  the  relationship  be¬ 
tween  the  standard  deviation  of  the  displacement  error  and  the 
mean  relative  modulus  estimation  error.  Fig.  9(b)  and  (d)  shows 
the  relationship  between  the  standard  deviation  of  the  displace¬ 
ment  error  and  the  standard  deviation  of  the  relative  modulus 
estimation  error. 
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(a)  (b) 


Fig.  9.  Modulus  estimation  performance  curves,  (a)  Standard  deviation  of  the  error  in  displacement  versus  the  mean  relative  error  in  modulus  estimates  resulting 
from  a  5%  standard  deviation  in  die  “measured”  force,  (b)  Standard  deviation  of  the  error  in  displacement  versus  standard  deviation  of  the  relative  error  in  modulus 
estimates  resulting  from  a  5%  standard  deviation  in  the  “measured”  force,  (c)  Standard  deviation  of  the  error  in  displacement  versus  the  mean  relative  error  in 
modulus  estimates  resulting  from  a  10%  standard  deviation  in  the  “measured”  force,  (d)  Standard  deviation  of  the  error  in  displacement  versus  standard  deviation 
of  the  relative  error  in  modulus  estimates  resulting  from  a  10%  standard  deviation  in  the  “measured”  force. 


Comparing  Fig.  9(a)  with  Fig.  9(c)  and  (d),  we  foimd  that  the 
modulus  estimation  error  is  not  veiy  sensitive  to  the  errors  in 
the  boundary  force  measurements.  However,  it  is  very  sensitive 
to  errors  in  the  displacement  measurements.  When  the  standard 
deviation  of  the  displacement  error  exceeds  10“^  mm,  the 
modulus  estimation  becomes  biased,  as  shown  in  Fig.  9(a) 
and  (b),  and  the  noise  in  modulus  estimation  starts  to  increase 
rapidly.  Note  that  the  units  in  the  horizontal  axis  of  all  plots 
in  Fig.  9  are  mm.  To  make  the  results  independent  of  the 
actual  object  dimension,  we  re-plotted  results  in  Fig.  10.  In 
Fig,  10  the  horizontal  axis  is  changed  to  standard  deviation  of 
relative  strain  error.  From  the  plots  shown  in  Fig.  10,  we  can 
see  that  when  the  noise  in  the  displacement  estimates  causes 
more  than  1%  strain  error,  the  quality  of  modulus  estimates 
starts  to  degrade  rapidly. 


IV,  Discussion 

The  derivation  in  Section  II-C2,  shows  that  the  boundary 
conditions  for  the  inverse  problem  with  the  proposed  method 
are  less  restrictive  than  those  of  the  forward  problem.  In 
the  forward  simulation,  boundary  conditions  that  describe  all 
sides  of  the  meshed  object  are  required  to  solve  (3)  for  the 
displacement  distribution  and  the  boundary  force  distribution. 
In  modulus  reconstruction,  however,  only  the  sub-ROI  of  the 
object  needs  to  be  meshed,  and  only  partial  boundary  force 
conditions  need  to  be  specified.  This  makes  our  approach  more 
practical  and  far  easier  to  implement  (experimentally)  than 
other  methods. 

In  our  method,  the  medium  is  assumed  to  be  elastic.  Possible 
tissue  viscous  behavior  is  not  accounted  for  in  our  model.  With 
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(a)  (b) 

Fig.  10.  Modiilus  estimation  performance  curves,  (a)  Standard  deviation  of  the  error  in  strain  versus  the  mean  relative  error  in  modulus  estimates  resulting  from 
a  5%  standard  deviation  in  Ae  ‘^measured”  force,  (b)  Standard  deviation  of  the  error  in  strain  versus  standard  deviation  of  the  relative  error  in  modulus  estimates 
resulting  from  a  5%  standard  deviation  in  die  “measured”  force,  (c)  Standard  deviation  of  the  error  in  strain  versus  the  mean  relative  error  in  modulus  estimates 
resulting  fi'om  a  10%  standard  deviation  in  the  “measured”  force,  (d)  Standard  deviation  of  the  error  in  strain  versus  standard  deviation  of  the  relative  error  in 
modulus  estimates  resulting  firom  a  10%  standard  deviation  in  the  “measured”  force. 


carefully  designed  methods  for  data  acquisition,  the  viscous  re¬ 
sponse  to  the  external  compression  can  be  negligible  [36]. 

In  Section  II-Cl  we  also  assume  the  elastic  behavior  of  the 
medium  to  be  linear.  This  assumption  only  needs  to  be  true  be¬ 
tween  two  states.  Si  and  52.  Most  human  tissues  have  a  non¬ 
linear  stress-strain  relationship.  Although  the  tissue  can  have 
significant  nonlinear  behavior,  we  can  restrict  our  analysis  to 
incremental  deformations  and  forces  as  described  in  (16).  Since 
the  incremental  deformation  between  Si  and  52  is  usually  small 
(less  than  2%),  the  assumption  of  linear  elasticity  is  reasonable. 

For  high  accuracy  FEA  solutions  to  the  forward  elasticity 
problem,  the  object  is  usually  meshed  with  nonuniform  size 
and  shape  elements.  The  mesh  has  higher  element  density 
near  curved  interfaces  where  the  modulus  changes  value,  but 


this  requires  knowledge  of  the  object  geometry.  In  modulus 
reconstruction,  the  internal  geometry  is  unknown,  so  rectan¬ 
gular  elements  are  used.  One  may  argue  that  it  is  possible 
to  first  obtain  a  rough  modulus  reconstruction  using  uniform 
elements,  then  re-mesh  the  object  with  nonvmiform  elements 
and  reconstruct  modulus  again  for  higher  accuracy.  To  do  so 
requires  higher  accuracy  displacement  estimates  in  the  regions 
with  higher  element  density.  Since  the  displacement  field  is 
usually  estimated  with  uniform  accuracy,  this  approach  is  likely 
to  fail. 

To  understand  the  effect  of  the  noise  in  the  measured  dis¬ 
placement  and  force  distribution,  we  conducted  a  number  of  nu¬ 
merical  simulations.  The  results  of  these  simulations,  shown  in 
Figs.  9  and  10,  provide  an  estimate  of  the  required  accuracy  in 
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input  displacement  and  force  data.  We  found  that  the  method  is 
less  sensitive  to  noise  in  the  force  measurements  than  to  noise 
in  the  displacement  estimates.  The  size  of  mesh  elements  can 
be  adjusted  to  change  the  sensitivity  to  displacement  estimation 
errors.  For  a  given  displacement  error,  larger  mesh  elements  re¬ 
sult  in  smaller  the  strain  error. 

From  Figs.  10(b)  and  (d),  we  observed  that  the  modulus  es¬ 
timates  obtained  from  noisy  displacement  and  force  estimates 
are  also  noisy  (22%  relative  errors).  However,  since  the  modulus 
contrast  between  normal  and  cancerous  tissues  is  usually  large 
(greater  than  100%)  [36],  the  modulus  image  contrast-to-noise 
ratio  for  cancerous  lesions  will  still  be  high. 

We  cannot  prove  that  the  inverse  problem  has  a  unique  solu¬ 
tion  from  physical  principles.  However,  our  method  provides  a 
unique  solution  algebraically.  For  the  simulations  that  we  have 
conducted,  when  there  is  no  noise  in  the  input  data,  the  solution 
that  we  obtained  is  the  same  (within  the  numerical  processing 
errors)  as  the  Young’s  modulus  distribution  that  we  specified  for 
the  object  (see  Fig.  7  and  Section  III-B).  With  the  added  noise 
(to  both  surface  force  and  the  displacement  distribution),  our 
method  generates  solutions  that  are  close  to  the  true  Young’s 
modulus  distribution(see  Fig.  9  and  Section  III-C).  This  sug¬ 
gests  that  our  method  is  stable  and  robust.  Since  we  lack  the 
necessary  equipment  to  simultaneously  measure  surface  force 
and  the  displacement  distribution  and,  therefore,  cannot  test  our 
method  experimentally.  This  will  be  the  subject  of  future  effort. 

FEA  treats  2-D  elasticity  problems  as  special  cases  of  a  gen¬ 
eral  three-dimensional  (3-D)  problem.  The  choice  of  these  spe¬ 
cial  cases  are  either  plain  strain  (elevational  strain  is  zero)  or 
plain  stress  (elevational  stress  is  zero).  With  the  plain  strain 
assumption,  the  external  force  is  assumed  to  be  exerted  on  a 
one-dimensional  boundary,  and  it  is  difficult  to  relate  such  a 
load  condition  to  reality.  With  the  plain  stress  assumption,  the 
object  has  finite  thickness,  and  the  calculated  boundary  force 
can  be  more  easily  related  to  actual  measurements  obtained  on 
a  2-D  surface. 

Tissue  deformation  is  3-D  in  nature.  However,  we  have  found 
that  in  vivo  breast,  for  example,  can  be  deformed  such  that  the 
motion  perpendicular  to  the  image  plane  is  small.  Thus,  a  2-D 
description  of  motion  provides  a  reasonable  approximation  to 
the  plane  strain  condition.  However,  force  measurements  are 
more  easily  related  to  the  plane  stress  condition,  and  with  this 
assumption  the  resulting  modulus  estimates  will  have  limited 
accuracy.  To  overcome  this  limitation,  we  need  to  extend  our 
approach  to  3-D.  Extending  2-D  modeling  to  3-D  is  relatively 
straightforward  with  FEA  methods. 

The  examples  and  discussion  of  displacement  estimation 
techniques  relate  to  our  work  in  using  ultrasound  to  track  tissue 
motion.  There  is  also  a  growing  body  of  work  where  magnetic 
resonance  techniques  are  used  for  estimating  tissue  elasticity 
[44],  [45].  The  modulus  reconstruction  technique  should  be 
applicable  in  that  work  as  well. 


reported.  Simulations  demonstrate  that  the  modulus  distribution 
for  an  ROI  can  be  determined  from  force  measurements  on  a 
single  surface  and  displacement  estimates  within  that  ROI.  The 
accuracy  in  force  and  displacement  estimates  required  with  this 
approach  are  also  estimated  with  simulations.  These  results  sug¬ 
gest  that  moduli  of  in  vivo  tissues  can  be  estimated  with  rea¬ 
sonable  accuracy  with  minor  modification  to  current  clinical 
imaging  systems. 


Appendix  I 

In  Fig.  1,  nodes  are  locally  numbered.^  The  element  nodal 
displacement  vector  for  the  rectangular  element  has  eight  com¬ 
ponents  and  can  be  written  as 


5(e) 


(4^ 


g(e) 


c(e)  c(e)  ^(e)  ^(e)  ^(e) 


^2x 


^2y 


^32/ 


(18) 


where  and  are  the  displacement  components  in  the  a: 
and  y  directions  of  the  first  node,  and  so  on  for  the  rest  of  the 
components;  T  is  the  matrix  or  vector  transpose  operator.  The 
element  nodal  force  vector  also  has  eight  components  and  can 
be  written  as 


Ae)  ,(e)  .(e)  .(e)  .(e)  .(e)  .(c) 

^  \J  lx  J  ly  /2x  J2y  J  Sy  J 4x  J4y  / 

(19) 

The  element  stiffhess  matrix  can  be  computed  as  [40] 

/b/2  .a/2  J, 

/  (20) 

-b/2  J-a/2 


where  t  is  the  element  thickness,  M  is  the  material  property  ma¬ 
trix,  and  B  is  the  displacement  to  strain  mapping  matrix.  For  ex¬ 
ample,  with  linear  elasticity  problems  and  assuming  plain  stress 
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where  and  are  the  Young’s  modulus  and  Possion’s 
ratio  of  element  e,  and  a  and  b  are  the  element  width  and  height. 
For  rectangular  elements,  the  displacement  to  strain  mapping 


matrix  is 
B^^'>  = 
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(22) 


V.  Conclusion 

A  new  approach  for  estimating  the  modulus  distribution  from 
noninvasively  determined  force  and  displacement  estimates  is 


^The  number  starts  from  the  lower  left  comer  of  the  element  and  increases  in 
the  clockwise  direction.  Note  that  other  numbering  methods  can  also  be  selected 
as  long  as  the  connectivity  matrix  (introduced  in  Step  2)  is  created  with  the  same 
numbering  scheme  for  all  elements. 
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where  i  =  1 , . , . ,  4  are  shape  functions.  The  shape  functions 
are  usually  selected  to  perform  bilinear  interpolations  that  have 
the  form 

The  representation  of  the  element  stiffness  matrix  (24)  can 
be  obtained  by  substituting  (21)-^23)  into  (20)  and  performing 
the  double  definite  integration.  For  convenience  in  deriving  the 
inverse  (modulus  estimation)  problem,  the  element  stiffness  ma¬ 
trix  can  be  written  as 

rM  =  (24) 

where  jFl  is  an  8-by-8  matrix  in  which  each  component  is  a  func¬ 
tion  of  the  aspect  ratio  (a/6)  and  Possion’s  ratio  of  the  rectan¬ 
gular  element. 

Appendix  II 

Following  is  a  description  of  the  global  stiffness  matrix  as¬ 
sembling  process. 

1)  Initialize  a  iV-by-iV  null  matrix  (all  zero  entries),  where 
N  equals  the  total  degrees  of  freedom  of  the  system  or 

N  =  number  of  nodes  x  degrees  of  freedom  per  node. 

(25) 

For  the  mesh  shown  in  Fig.  2,  A  =  16  x  2  =  32. 

2)  For  element  e,  generate  a  local  (element)  variable  number 
to  global  (system)  variable  number  conversion  index 
vector 

)  =  ^Cei*2  1,  Cei+2,  Ce2*2  1,  Ce2*2,  Ce3’i'2  1, 

Ue3*2,  Ce4*2  1,  Ce4*2^  (26) 

where  Cei  is  the  row  and  first  column  entry  of  the 
connectivity  matrix  C,  and  so  on.  For  the  system  shown 
in  Fig.  2,  the  index  vector  for  the  first  element  is 

=  2  3  4  11  12  9  10).  (27) 

3)  Accumulate  to  Kjie)jie)  for  i  =  1, . . . ,  8  and  j  = 

1, . . . ,  8  (K^f  is  the  ith  row  jth  column  component  of 
matrix  is  the  row  column  com¬ 
ponent  of  matrix  K ;  is  the  «th  component  of  the  index 

matrix). 

4)  Iterate  2  and  3  for  all  elements. 

Appendix  HI 

The  displacement  boundary  condition  can  be  defined  as 
Ssuh  ~  where  6sub  is  a  vector  that  is  composed  of  a  subset 


of  the  components  of  the  global  nodal  displacement  vector  6; 
Sc  is  a  known  constant  vector  that  specifies  the  nodal  boundary 
displacement.  For  example,  a  common  displacement  boundary 
condition  for  the  object  meshed  by  our  9-node  example  shown 
in  Fig.  2  is  to  compress  the  top  side  downward  1%  of  the 
total  height  of  the  object  while  the  bottom  is  fixed  vertically. 
This  example  displacement  boundary  condition  can  be  defined 
as  ^sub  =  (^2  ^4  ^6  ^8  ^26  ^28  ^30  ^32)^  and 

Sc  =  (0  0  0  0  O.Ol/i  O.Ol/i  O.Ol/i  0.01ft  )^, 

where  h  is  the  height  of  the  object.  The  penalty  approach  can 
be  expressed  as  the  following  seven  steps. 

1)  Initialize  the  global  force  vector  /  as  a  null  vector. 

2)  Select  a  large  number  L  (a  choice  for  L  is  L  =  max\K\  x 
10^  as  suggested  by  Chandmpatla  [41], 

3)  According  to  6sub,  set  the  corresponding  /sub  to  L  x  6c 

(/sub  =  ( /2  /4  /e  /s  /26  /28  ho  fz2  for 

the  given  example). 

4)  According  to  6s  ub,  add  L  to  the  corresponding  diagonal 

component  of  K  (for  our  example,  these  diagonal  compo¬ 
nents  are  Ff2,2,  -^4,4,  ^6,6,  26)  -^^28,28)  ^30,30) 

and  ^^32, 32)- 

5)  Solve  iif6  =  /  for  6. 

6)  Calculate  reaction  force  /sub  =  —L{Ssuh  —  Sc). 

7)  Replace  6sub  with  6c. 

References 

[1]  C.  P.  McPherson,  K,  K.  Swenson,  G.  Jolitz,  and  C.  L.  Murray,  “Survival 
among  women  ages  40-49  years  with  breast  carcinoma  according  to 
method  of  detection,”  Cancer,  vol.  79,  no.  10,  pp.  1923-1932,  1997. 

[2]  A.  Sarvazyan,  “Mechanical  imaging:  A  new  techonology  for  medical 
diagonsis,”  Int  J.  Med.  Inform.,  vol.  49,  pp.  195-216,  1998. 

[3]  P.  S,  Wellman,  ‘Tactile  imaging,”  Ph.D.  dissertation.  Dept,  Eng.  Sci., 
Harvard  Univ.,  Cambridge,  MA,  1999. 

[4]  E.  I.  Cespedes,  J.  Ophir,  H,  Ponnekanti,  and  N.  Maklad,  “Elastography: 
Elasticity  imaging  using  ultrasound  with  application  to  muscle  and 
breast  in  vivof  Ultrasound  Imag.,  vol.  15,  pp.  73-88, 1993. 

[5]  E,  Chen,  R.  Adler,  P,  Carson,  W.  Jenkins,  and  W.  O’Brein,  “Ultra¬ 
sound  tissue  displacement  imaging  using  ultrasound  with  application 
to  muscle  and  breast  cancer,”  Ultrasound  Med.  Biol.,  vol.  21,  pp. 
1153-1162,  1995. 

[6]  B.  S.  Garra,  E.  1.  Cespedes,  J.  Ophir,  S.  R.  Spratt,  R.  A.  Zuubrier,  C.  M. 
Magnant,  and  M.  F.  Pennanen,  “Elastography  of  breast  lesions:  Initial 
clinical  results,”  Radiology,  vol.  202,  pp.  79-86,  1997, 

[7]  T.  J.  Hall,  Y.  Zhu,  and  C.  S.  Splading,  “/w  vivo  real-time  freehand  palpa¬ 
tion  imaging,”  Ultrasound  Med.  Biol,  Oct  2002. 

[8]  S.  Y.  Emelianov,  M.  A.  Lubinski,  W.  F.  Weiztel,  R.  C.  Wiggins,  A.  R 
Skovoroda,  and  M.  O’Donnell,  “Elasticity  imaging  for  early  detection  of 
renal  pathology,”  Ultrasound  Med.  Biol,  vol.  21,  pp.  871-883, 1995. 

[9]  P,  Chatruvedi,  M.  F.  Insana,  and  T.  J.  Hall,  “Acoustic  and  elastic  imaging 
to  model  disease-induced  changes  in  soft  tissue  structure,”  in  Lecture 
Notes  in  Computer  Science,  A,  Kuba  and  M.  Samal,  Eds,  New  York: 
Springer-Verlag,  1997,  vol.  1230,  pp.  1-14. 

[10]  F.  Kallel,  J.  Ophir,  K.  McGee,  and  T.  Krouskop,  “Elastographic  imaging 
of  low-contrast  elastic  modulus  distributions  in  tissue,”  Ultrasound  Med. 
Biol,  vol.  24,  pp.  409-425, 1998, 

[1 1]  T.  J.  Hall,  H.  A.  Khant,  M.  F.  Insana,  J.  G.  Wood,  Y.  Zhu,  D.  Preston, 
and  B.  D,  Cowley,  “The  utility  of  quantitative  ultrasound  for  tracking 
the  progression  of  polycystic  kidney  disease,”  Proc.  SPIE,  vol.  3982, 
pp.  161-171, 2000. 

[12]  B.  M.  Shapo,  J.  R.  Crowe,  A.  R,  Skovoroda,  M.  J.  Eberle,  N.  A.  Cohn, 
and  M.  O’Donnel,  “Displacement  and  strain  imaging  of  coronary  ar¬ 
teries  with  intraluminal  ultrasound,”  IEEE  Trans.  Ultrason.,  Ferroelecl 
Freq.  Contr.,  vol.  43,  pp.  234-246,  Mar.  1996. 

[13]  L.  K.  Ryan  and  F,  S.  Foster,  “Ultrasonic  measurement  of  differential 
displacement  and  strain  in  a  vascular  model,”  Ultrason.  Imag.,  vol.  19, 
pp.  19-38,  1997. 

[14]  Y.  C.  Fung,  A  First  Course  in  Continuum  Mechanics.  Englewood 
Cliffs,  NJ:  Prentice-Hall,  1994. 


ZHU  et  al.:  A  FINITE-ELEMENT  APPROACH  FOR  YOUNG’S  MODULUS  RECONSTRUCTION 


901 


[15]  T-  A.  Krouskop,  D.  R.  Dougherty,  and  F.  S.  Vinson,  “A  pulsed  doppler 
ultrasonic  system  for  making  noninasive  measurement  of  the  mechanical 
properties  of  soft  tissue,”  J.  Rehab.  Res.  Dev.,  vol.  24,  pp.  1-8, 1987. 

[16]  R.  M.  Lemer,  S.  R.  Huang,  and  K.  J.  Parker,  “Sonoelasticity  images  de¬ 
rived  from  ultrasound  signals  in  mechanically  vibrated  tissues,”  Ultra¬ 
sound  Med.  Bioi,  vol.  16,  no.  3,  pp.  231-239, 1990. 

[17]  Y.  Yamakoshi,  J.  Sato,  and  T.  Sato,  “Ultrasonic  imaging  of  internal  vi¬ 
bration  of  soft  tissue  under  forced  vibration,”  IEEE  Trans.  Ultrason., 
Ferroelect.  Freq.  Contr.,  vol.  37,  pp.  45-53,  Mar  1990. 

[18]  L.  Sandrin,  S.  Catheline,  M.  Tanter,  X.  Hannequin,  and  M.  Fink, 
“Time-resolved  pulsed  elastography  with  ultrafast  ultrasonic  imaging,” 
Ultrason.  Imag.,  vol.  21,  no.  4,  pp.  259-272,  1999. 

[19]  R.  J.  Dickinson  and  C.  R.  Hill,  “Measurement  of  soft  tissue  mation  using 
correlation  between  a-scans,”  Ultrasound  Med.  Biol,  vol.  8,  no.  3,  pp. 
263-271, 1982. 

[20]  L.  S,  Wilson  and  D.  E.  Robinson,  “Ultrasonic  measurement  of  small 
displacements  and  deformations  of  tissue,”  Ultrason.  Imag.,  vol.  4,  pp. 
71-82,  1982. 

[21]  M.  Tristam,  D.  C.  Barbosa,  D.  O.  Cosgrove,  D.  K.  Nassiri,  J.  C.  Bamber, 
and  C.  R.  Hill,  “Ultrasonic  study  of  in  vivo  kinetic  characteristics  of 
human  tissues,”  Ultrasound  Med.  Biol,  vol.  12,  no.  12,  pp.  927-937, 
1986. 

[22]  J.  Ophir,  1.  Cespedes,  H.  Poimekanti,  Y  Yazdi,  and  X.  Li,  “Elastography: 
A  quantitive  method  for  imaging  the  elasticity  of  biological  tissues,” 
Ultrason.  Imag.,  vol.  13,  pp.  111-134,  1991. 

[23]  M.  O’Donnell,  A.  Skovoro^  B.  Shapo,  and  S.  Emelianov,  “Internal  dis¬ 
placement  and  strain  imaging  using  ultrasonic  speckle  tracking,”  IEEE 
Trans.  Ultrason.,  Ferroelect.  Freq.  Contr.,  vol.  41,  pp.  314-325,  May 
1994. 

[24]  M,  Bilgen  and  M.  F,  Insana,  “Deformation  models  and  correlation  anal¬ 
ysis  in  elastography,”  J.  Acoust.  Soc.  Amer.,  vol.  99,  pp.  3212-3224, 
1996. 

[25]  K.  A.  Wear  and  R.  L.  Popp,  ‘Theorectical  analysis  of  a  technique  for 
die  characterization  of  myocardium  contraction  based  upon  temporal 
correlation  of  ultrasound  echoes,”  IEEE  Trans.  Ultrason.,  Ferroelect. 
Freq.  Contr.,  vol.  UFFC~34,  pp.  368-375,  1987. 

[26]  R.  Adler,  J.  Rubin,  P.  Bland,  and  P.  Carson,  “Quantitive  tissue  motion 
analysis  of  digitized  m-mode  images:  Gestational  differences  of  fetal 
lung,”  Ultrasound  Med.  Biol,  vol.  16,  pp.  561-569,  1990. 

[27]  P.  DeJong,  T.  Arts,  A.  Hoeks,  and  R.  Reneman,  “Determination  of  tissue 
motion  velocity  by  correlation  interpolation  of  pulsed  ultrasonic  echo 
signals,”  Ultrason.  Imag.,  vol.  12,  pp.  84-98, 1990, 

[28]  Y  Zhu  and  T.  J,  Hall,  “A  modified  block  matching  metiiod  for  realtime 
freehand  strain  imaging,”  Ultrason.  Imag.,  Nov.  2002. 

[29]  R-  Muthupillai  and  R.  L.  Ehman,  “Magnetic  resonance  elastography,” 
Nat.  Med.,  vol.  2,  no.  5,  pp.  601-603,  1996. 

[30]  A.  M.  T.  E.  Oliphant,  M.  A.  Dresner,  J.  L.  Mahowald,  S.  A.  Kruse,  E. 
Axnromin,  J.  P.  Fehnlee,  J.  F.  Greenleaf,  and  R.  L.  Ehman,  “Magnetic 
resonance  elastography;  Non-invasive  mapping  of  tissue  elasticity,” 
Med.  Image  Anal,  vol.  5,  no.  4,  pp.  237-254, 2001. 


[31]  J.  B.  Weaver,  E.  E.  V.  Houten,  M.  I.  Miga,  F.  E.  Kennedy,  and  K  D. 
Paulsen,  “Magnetic  resonance  elastography  using  3D  gindient  echo 
measurements  of  steady-state  motion,”  Med.  Phys.,  vol,  28,  no,  8,  pp. 
1620-1628,2001. 

[32]  D.  D.  Duncan  and  S .  J.  Kirkpartrick,  ‘Trocessing  algorithms  for  tracking 
speckle  shifts  in  optical  elastography  of  biological  tissues,”,/  Biomed. 
Opt.,  vol.  6,  no.  4,  pp.  418-426,  2001. 

[33]  A.  R.  Skovoroda,  S.  Y.  Emelianov,  and  M.  O’Donnell,  “Tissue  elasticity 
reconstruction  based  on  ultrasonic  displacement  and  strain  images,” 
IEEE  Trans.  Ultrason.,  Ferroelect.  Freq.  Contr.,  vol.  42,  pp.  747-765, 
July  1995. 

[34]  F.  Kallel  and  M.  Bertrand,  “Tissue  elasticity  reconstruction  using  linear 
perturbation  method,”  IEEE  Trans.  Med.  Imag.,  vol.  15,  pp.  299-313, 
June  1996. 

[35]  M.  M,  Doyley,  P.  M.  Meaney,  and  J.  C.  Bamber,  “Evaluation  of  an  iter¬ 
ative  reconstruction  method  for  quantitative  elastography,”  Phys.  Med. 
Biol,  vol.  45,  pp.  1521-1540, 2000. 

[36]  T.  A.  Krouskop,  T.  M.  Wheeler,  F.  Kallerl,  B.  S.  Garra,  and  T.  J.  Hall, 
“Elastic  moduli  of  breast  and  prostate  tissues  under  compression,”  Ul¬ 
trason.  Imag.,  vol.  20,  pp.  260-274, 1998. 

[37]  K.  R-  Raghavan  and  A.  E.  Yagle,  “Forward  and  inverse  problems  in 
elasticity  imaging  of  soft-tissue,”  IEEE  Trans.  Nucl  Set,  vol.  41,  pp. 
1639-1648,  Aug.  1994. 

[38]  A.  J.  Romano,  J.  J.  Shirron,  and  J.  A.  Bucaro,  “On  the  noninvasive  deter¬ 
mination  of  material  parameters  from  a  knowledge  of  elastic  displace¬ 
ments:  Theory  and  numerical  simulation,”  IEEE  Trans.  Ultrason.,  Fer¬ 
roelect.  Freq.  Contr.,  vol.  45,  pp.  751-759,  May  1998. 

[39]  K.  H.  Huebner,  E.  A.  Thornton,  and  T.  G.  Byrom,  The  Finite  Element 
Method  for  Engineers,  3rd  ed.  New  York:  Wiley,  1995. 

[40]  Y  W.  Kwon  and  H.  Bang,  The  Finite  Eement  Method  Using  MATLAB, 
2nded.  New  York:  CRC,  2000. 

[41]  T.  R.  Chandrupatla  and  A.  D.  Belegundu,  Introduction  to  Finite  Elements 
in  Engineering,  2nd  ed.  Englewood  Cliffs,  NJ:  Prentice-Hall,  1997. 

[42]  W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flaimery, 
“Numerical  Recipes  in  C,”  in  The  Art  of  Scientific  Computing,  2nd 
ed.  Cambridge,  MA:  Cambridge  Univ.  Press,  1992. 

[43]  M.  Bilgen  and  M.  F.  Insana,  “Convariance  analysis  of  time  delay  esti¬ 
mates  for  strained  signals,”  IEEE  Trans.  Signal  Processing,  vol.  46,  pp. 
2589-2600,  Oct.  1998. 

[44]  E.  E.  V.  Houen,  K.  D.  Paulsen,  M.  I.  Miga,  F.  E.  Kennedy,  and  J.  B. 
Weaver,  “An  overlapping  subzone  technique  for  MR-based  elastic  prop¬ 
erty  reconstruction,”  Magn.  Reson.  Med.,  vol.  42,  no.  4,  pp.  779-786, 
1999. 

[45]  T.  Wu,  J.  P,  Felmlee,  J.  F.  Greenleaf,  S.  J.  Riederer,  and  R.  L.  Ehman, 
“Assessment  of  thermal  tissue  ablation  with  MR  elastography,”  Magn. 
Reson.  Med.,  vol.  45,  no.  1,  pp.  80-87,  2001. 


Noise  Reduction  Strategies  in  Freehand  elasticity  Imaging 


Timothy  J.  Hall,  Jingfeng  Jiang,  Yanning  Zhu,  Lany  T  Cook 
Department  of  Radiology,  University  of  Kansas  Medical  Center 
3901  Rainbow  Boulevard,  Kansas  City,  KS  66160-7234 


Abstract —  We  are  developing  a  clinical  ultrasonic 
imaging  system  for  real-time  estimation  and  display 
of  tissue  elastic  properties.  We  have  demonstrated 
that  real-time  feedback  of  elasticity  images  is  essen¬ 
tial  for  obtaining  high-quality  data  (consecutive  im¬ 
ages  with  high  spatial  coherence).  The  key  element  to 
successful  scanning  is  real-time  visual  feedback  y/hich 
guides  die  patient  positioning  and  compression  direc¬ 
tion.  Our  data  have  clearly  demonstrated  nonlinear¬ 
ity  in  the  strain  properties  of  different  tissue  types. 
We  have  also  demonstrated  that  a  comparison  of  the 
area  of  a  breast  lesion  observed  in  strain  images  ver¬ 
sus  B-mode  images  is  a  sensitive  criterion  for  differ¬ 
entiating  malignant  from  benign  tumors.  Frame-to- 
frame  variability  in  strain  images  somewhat  degrades 
the  ability  to  observe  diese  phenomena.  Three  strate¬ 
gies  for  reducing  frame-to-fi^e  strain  image  noise  are 
described.  The  combination  of  these  post-processing 
strate^es  provides  a  signicant  improvement  in  the 
quality  of  long  sequences  of  strain  images. 

I.  Introduction 

We  are  implementing  and  testing  real-time  mechan¬ 
ical  strain  imaging  integrated  into  a  clinical  ultrasound 
imaging  system  (Elegra,  Siemens  Medical  Solutions) 
{!].  Our  work  was  motivated  by  promising  in  vivo  re¬ 
sults  reported  by  Garra  et  al.  [2].  In  that  report  they 
described  data  acquisition  based  on  a  modied  mam¬ 
mography  system-  The  use  of  that  system  limited  the 
areas  of  the  breast  from  which  they  could  acquire  data. 

We  are  grateful  for  the  nancial  support  of  USAMRAA 
DAMDl  7-00-1-0596  and  technical  support  from  Siemens  Med¬ 
ical  Systems  Ultrasound  Group.  The  U.S.  Army  Medical  Re¬ 
search  Acquisition  Activity,  820  Chandler  Street,  Fort  Detrick  MD 
21702-5014  is  the  awarding  and  administering  acquisition  ofce 
for  DAMD 17-00-1 -05%.  The  information  reported  here  does  not 
necessarily  rcect  the  position  or  policy  of  the  U.S.  Government, 
and  no  ofeial  endorsement  should  be  inferred. 
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Freehand  scanning  has  been  the  dominant  method  of 
clinical  sonography  for  many  years.  So,  freehand  scan¬ 
ning  will  likely  more  quickly  gain  clinical  acceptance 
of  elasticity  imaging  if  it  can  be  performed  efciently  , 
We  have  argued  that  real-time  feedback  to  the  hand- 
eye  coordination  system  allows  constant  manipulation 
of  the  boundary  conditions  of  deformation  and  allows 
die  observer  to  know  when  high  quality  strain  image 
data  are  acquired.  The  small  delay  between  acquiring 
successive  frames  (tens  of  milliseconds)  and  the  rela¬ 
tively  slow  deformation  rate  (cyclic  freehand  deforma¬ 
tion  at  about  1  Hz)  likely  results  in  primarily  an  elastic 
response  in  tissue  (minimal  viscous  effect)  [3]. 

In  vivo  elasticity  images  of  breast  lesions  obtained 
with  our  system  have  high  contrast-to-noisc  ratios,  hi 
fact,  relatively  long  sequences  (30  sequential  frames 
or  more)  of  high  quality  strain  images  are  normally 
obtained  in  clinical  trials.  An  example  image  from  a 
breast  tumor  is  shown  in  Fig.  1.  Using  our  system 
we  have  demonstrated  that  broadenomas  often  have 
a  surface  pressure-dependent  strain  image  contrast  [IJ. 
The  strain  image  contrast  for  a  broadenoma  is  gener¬ 
ally  highest  with  the  least  surface  pressure  and  contrast 
(tecreases  as  the  pressure  is  increased.  In  addition,  the 
relative  size  of  a  lesion  in  B-mode  versus  strain  im¬ 
ages  is  a  sensitive  criterion  for  differentiating  malig¬ 
nant  from  benign  lesions  [1,2].  A  plot  of  lesion  area 
measured  in  B-mode  Images  versus  lesion  area  in  the 
corresponding  strain  images  (Fig.  2)  illustrates  this  per¬ 
formance. 

Viewing  strain  image  sequences  allows  visualiza¬ 
tion  of  strain  image  noise  that  is  not  apparent  in  single 
strain  images.  Given  that  pressure  dependent  strain  im¬ 
age  contrast  and  the  relative  lesion  size  comparison  are 
useful  criteria  for  (hfferentiating  breast  lesions,  reduc¬ 
ing  that  frame-to-frame  variability  in  strain  estimates 
will  improve  the  ability  to  visualize  that  vaiying  con¬ 
trast  and  better  determine  lesion  boundaries. 
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Fig.  1.  A  B'mode  strain  image  pair  obtained  by  freehand 
Bcanning  of  a  breast  cardnoma  in  vivo. 


Fig.  X  A  companscm  of  lesion  area  for  carcinomas  (circles),  cysts 
(squares)  and  fibroadenomas  (diamonds). 


Stxat^es  for  xeducmg  that  &anie-to<-fraixie  variabil¬ 
ity  in  local  ^rain  estimates  is  the  subject  of  this  report. 
The  goal  is  to  reduce  noise  while  tnaintainiag  spatial 
resolution.  An  additional  goal  is  to  obtain  reduction  in 
strain  image  noise  with  minimal  incremental  computa- 
donai  load  such  that  the  technique  can  be  implemented 
in  (Ht-line  [SDcessiiig  for  clinical  trials.  Several  pos¬ 
sible  approaches  are  suggested,  and  results  of  imple¬ 
menting  those  approaches  are  encouraging. 

IL  Materials  and  Methods 
Strain  Image  Formation 

A  2-D  block  matching  algorithm,  based  on  the  sum- 
squared  difference  (SSD)  algorithm,  is  used  for  mo¬ 
tion  tracking  in  our  implemencadon.  The  kernel  size 
was  selected  to  approximate  the  2-D  pulse-echo  ui^ 
trasound  point  Sfnead  function  for  the  system  em¬ 
ployed  (Siemens  SONOLINE  Hlegra  with  75L40  and 
VFX13-5  linear  arrays).  Data  were  processed  on  the 


image  processor  subsystem  of  the  Elegra.  The  algo¬ 
rithm  displays  streaming  B-mode  and  strain  images 
side-by-side  at  about  seven  frames  per  second  and 
stares  the  full  sequence  of  I-Q  echo  data  for  on-line 
post-processing. 

Patient  Scawung 

All  patients  provided  informed  consent  consistent 
with  the  protocol  approved  by  the  Human  Subjects 
Committee  (Institutional  Review  Board)  at  Kansas 
University  Medical  Center.  Patient  scans  were  per- 
foimed  in  a  manner  consistent  with  a  normal  breast 
ultrasound  exam;  the  breast  was  scanned  with  the  pa¬ 
tient  (typically)  in  the  supine  position  with  her  ann  be¬ 
hind  her  head.  When  the  Iweast  lesion  was  located, 
the  transducer  was  pressed  toward  the  chest  wall  at 
a  steady  rate  in  an  effort  to  achieve  about  0.5-1.2% 
compression  frame-to-frame  while  repeating  the  com- 
piess/release  cycle  for  relatively  large  (>10%)  com¬ 
pression.  The  compression  technique  was  adjusted,  by 
changing  the  compression  direction  or  patient  position, 
until  there  was  nearly  uniaxial  compression  with  min¬ 
imal  elevation  motion.  Real-time  B-mode  and  strain 
image  display  allowed  visualization  of  the  data  qual¬ 
ity.  Using  this  scanning  technique,  no  patient  has  ex¬ 
perienced  any  discomfort  in  our  procedures. 

III.  Sources  of  Noise  and  Strategies  for 
Improved  results 

Small  Deformation-Law  Average  Strain 

The  strain  image  contrast-to-noise  ratio  increases 
as  the  applied  (uniaxial)  strain  increases  for  fcarae- 
av€^ge  strains  up  to  about  5%  (see,  for  example  [4]), 
Freehand  deformation  sometimes  results  in  nonuni- 
form  firame-to-frame  average  strain.  The  pairing  of  if 
echo  frames  of  dau  in  post-processing  can  be  adjusted 
to  select  appropriate  frame  pairs  to  obtain,  nominally, 
1-1.5%  frame-average  strain  [5].  A  representative  re¬ 
sult  for  a  m  vivo  breast  fibroadenoma  is  shown  in  fig¬ 
ure  3. 

Large  Displacement  Errors 

Regardless  of  the  average  applied  strain,  the  motion 
tracking  algorithm,  and  whether  that  strain  was  applied 
with  freehand  scanning  or  motorized  deformation,  dis¬ 
placement  estimation  errors  sometime  occur.  When 
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Fig.  3.  The  fnime'lD>firame  strain  obtained  by  freehand  scan^ 
Ding  of  a  breast  ir  vivo  and  pairing  adjacent  frames  (dashed)  or 
dynamically  adjusting  the  frame  pidriog  in  post-processing  (solid) 
for  stridn  estsmalkm. 


Fig.  4.  Displacement  estimates  resulting  from  2-D  block  matching 
(dast^)  and  moving  linear  regression  (solid)  for  those  results  as  a 
Unction  of  d^th  for  one  line  of  data. 


using  a  very  small  2-D  kernel  to  trade  motion,  rela¬ 
tively  large  enors  are  easily  detected  A  statistical  ar¬ 
gument  and  a  moving  linear  regression  can  be  used  to 
detect  displacement  estimates  that  ate  so  different  from 
their  neighbors  that  they  are  very  unlikely  to  be  cot- 
lea.  Foe  exanqile,  if  we  know  thsx  adjacent  displace¬ 
ment  estimates  diff^ences  of  one  sample  imply  a  local 
strain  of  7%,  and  the  frame-average  strain  is  only  about 
1%,  then  adjacent  displacement  estimate  differences  of 
more  than  one  san^le  are  extremely  uhlilcely.  We  can 
use  lineariegtession  to  estimate  the  local  displacement 
and  ccHi^are  each  ei^imate  to  the  local  regression  fit.  If 
tbt  diffimnce  between  the  individual  estimate  and  the 
regression  value  is  more  than  some  threshold,  ex- 
ample  one  sample,  that  estimate  is  judged  to  be  in  error 
and  is  replaced  by  the  regression  \^ue.  This  approach 
is  very  ^ective  at  detecting  and  removing  large  dis¬ 
placement  mors  prkHT  to  estiinattiig  local  strain.  An 
example  of  this  technique  applied  to  in  vivo  breast  data 
is  shown  in  Fig.  4 

It  is  straightforward  to  extend  this  approach  and  fit 
local  displaconem  estimates  to  a  small  (planar)  sur¬ 
face  to  d^ect  and  exdude  large  errors.  The  advantage 
of  using  a  surface  fit  is  that  it  minimizes  local  axial 
shean  The  size  of  the  surface,  as  with  the  linear  re¬ 
gression  window,  is  a  tradeoff  between  tracking  true 
di^acements  (resolution)  with  small  surfaces  and  re¬ 
jecting  consecutive  displacement  errors  (noise  reduc¬ 
tion)  with  a  longer  window,  A  2  mm  x  2  mm  surface 
provides  subjectively  satisfactory  results. 


Complex  Motion 

Correlation  search  techniques  for  motion  tracking, 
such  as  that  employed  for  our  real-time  strain  imag¬ 
ing,  woric  well  when  the  deformation  is  small  and 
the  motion  is  relatively  sitrple  [6].  When  the  mo¬ 
tion  is  relatively  large  (>5%  strain)  or  includes  obvi¬ 
ous  shear  the  performance  of  correlation  search  tech¬ 
niques  degrades.  Under  these  conditions  higher  order 
motion  models  ^ically  improve  motion  tracking  per¬ 
formance  but  at  higher  computational  load  [7]. 

In  these  techniques,  a  cost  function  is  defined  that 
minimizes  the  constramts  used  to  model  the  motion. 
The  basic  approadi  assumes  that  the  brightness  of  each 
point  in  the  image  remains  constant  from  one  frame  of 
data  to  the  next  A  brightness  constraint,  is  then 
used  to  estimate  the  defotmation  field  that  describes 
the  transformation  from  one  firame  of  data  to  the  next 
That  general  concept  was  adapted  to  rf  echo  samples, 
Instead  of  B-mode  pixd  brightness  in  our  previous 
work  [8].  Our  current  approach  uses  the  sum  squared 
dififeience  between  pie-  and  warped  post-compressiOD 
echo  fields  as  the  brightness  constraint  A  smoothness 
constraint,  is  added  to  the  cost  function  to  place 
a  penalty  on  large  deviations  from  average  local  mo¬ 
tion.  Numerous  approaches  to  create  smoothness  con¬ 
straints  have  been  reported.  One  of  the  earliest  smooth¬ 
ness  constraints  was  reported  by  Horn  [9]. 

From  our  previous  work  using  similar  techniques  [8] 
we  recognized  the  importance  of  the  constraint  func¬ 
tion  in  the  overall  performance  of  this  approach.  We 
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ate  investigating  the  use  of  strain  eneigy  as  a  smooth¬ 
ness  constraint  function,  as  suggested  by  Bajscy  [lO]. 
The  resulting  cost  function  is 

C  =  jj{yEi  +  E,)dxdy  (1) 

a 

where  7  is  an  adaptively  chosen  scale  factor.  Example 
images  of  the  same  data  with  and  without  these  noise 
reduction  strategies  are  shown  in  Figs.  5  and  6. 


Fig.  5.  A  B  mode  and  strain  image  pair  without  noise  reduction. 


Fig.  6.  A  B-mode  and  strain  im^c  pair  with  all  three  noise  re¬ 
duction  strat^es. 


IV.  Discussion 

Real-time  display  of  side-by-side  B-mode  and  strain 
images  is  essential  for  guiding  the  manipulation  of 
boundary  conditions  for  the  mechanics  experiment  that 
is  strain  imaging.  The  real-time  feedback  to  the  hand- 
eye  coordination  systems  allows  the  sonographer  to 
manipulate  the  compression  direction,  force,  and  rate 
to  obtain  high-quality  sequences  of  strain  images.  The 
system  involves  no  remote  data  acquisition  or  display, 
and  no  additional  ^gnal  processing  hardware.  It  is 
fully  integrated  into  the  Elegra  system. 

The  sequence  of  B-mode  and  strain  image  pairs  al¬ 
lows  the  sonographer  to  select  images  representative 


of  the  “typical”  strain  image  for  a  lesion.  This  abil¬ 
ity,  along  with  better  determination  of  lesion  bomkl- 
ary  available  by  viewing  a  sequence  of  images,  has 
likely  improved  the  ability  to  measure  true  lesion  size 
in  strain  imaging  compared  with  the  results  reported  by 
Garra,  et  al. 

Deformable  model  approaches  are  iterative  and 
require  an  initial  displacement  field  approximation. 
When  that  approximation  has  large  errors,  accurate 
displacement  estimates  axe  influenced  by  the  presence 
of  local  displacement  errors.  By  first  using  the  above 
techniques  to  minimize  local  displacement  errors,  we 
start  the  deformable  model  iteration  with  a  reasonable 
approximation  to  the  true  displacement  field.  The  re¬ 
sult  is  a  significam  reduction  in  firamc-to-frame  vari¬ 
ability  in  strain  images  and  improved  detection  of  le¬ 
sion  boundaries  and  contrast  changes. 

V.  Conclusions 

Strategies  for  reducing  strain  image  noise  in  a  se¬ 
quence  of  strain  images  are  described.  Tests  of  these 
strategies  on  breast  data  acquired  in  vivo  demonstrate 
improvements  in  the  image  sequence.  In  combination 
these  strategies  provides  significant  improvements  in 
the  visual  interpretation  of  a  sequence  of  strain  images. 
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Abstract 


A  new  mode  of  imaging  with  ultrasoimd  is  under  development  in  several  laboratories 
around  the  world.  The  new  mode  estimates  some  measure  of  the  viscoelastic  properties  of 
tissue.  The  information  displayed  in  these  images  is  a  surrogate  for  that  obtained  with 
manual  palpation.  This  report  reviews  some  of  the  fuiidamental  concepts  in  elasticity 
imaging  and  highlights  the  development  of  a  system  for  a  real-time  elasticity  imaging 
being  tested  for  in  vivo  breast  and  thyroid  imaging.  Results  of  early  laboratory  tests  that 
motivated  this  effort  are  reviewed  and  results  obtained  using  this  system  with  in  vivo 
tissues  are  included.  Imaging  the  elasticity  of  in  vivo  breasts  suggest  that  invasive  ductal 
carcinomas  appear,  on  average,  more  than  twice  as  large  in  an  elasticity  image  compared 
to  the  same  lesion  in  an  ultrasound  B-mode  image,  but  fibroadenomas  and  cysts  are 
nearly  equal  in  size  in  the  two  image  types.  Data  from  that  real-time  system  also 
demonstrate  that  the  relative  stiffness  of  many  fibroadenomas  change  as  they,  and  their 
surrounding  tissue,  are  deformed.  The  utility  of  this  technology,  and  the  new  information 
it  provides,  suggests  that  it  might  soon  be  available  on  commercial  ultrasoimd  imaging 
systems. 


Summary  Statement 

Elasticity  imaging  could  significantly  improve  the  diagnosis  of  breast  abnormalities  using 
ultrasound. 
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INTRODUCTION 

It  is  the  experience  of  many  that  palpation,  pressing  on  the  surface  of  soft  tissue  in 
an  effort  to  ‘feel’  abnormalities,  is  a  commonly  used  diagnostic  tool.  This  tool  has  been 
used  for  thousands  of  years  and  is  the  primary  diagnostic  tool  for  some  diseases. 
Examples  include  breast  self  examination  for  sensing  breast  ‘lumps’  and  digital  rectal 
examination  for  prostate  cancer.  Palpation  is  known  to  be  subjective  and  it  lacks 
sensitivity  to  small  abnormalities  that  are  deep  beneath  the  skin  surface. 

Improving  sensitivity  and  reducing  the  subjectivity  of  palpation  could  have  a 
significant  impact  on  breast  cancer  prognosis.  Breast  cancer  is  the  second-leading  cause 
of  cancer  deaftis  in  women.  Over  200,000  new  cases  of  invasive  breast  cancer  are 
expected  in  the  USA  this  year  alone.  It  is  anticipated  that  approximately  40,000  women 
in  the  USA  will  die  of  breast  cancer  in  2003.  The  prognosis  for  breast  cancer  patients  is 
best  when  the  disease  is  detected  at  an  early  stage.  Specifically,  5-  and  10-year  survival 
statistics  are  best  when  cancer  is  noninvasive’  and  is  less  than  1cm  in  diameter.^ 
Improvements  in  mammography  have  resulted  in  iniproved  detection  of  breast  lesions, 
and  mammography  has  been  shown  capable  of  detecting  smaller  tumors  in  yoimg  women 
than  either  breast  self-examination  or  clinical  breast  examination.^  However, 
mammography  is  not  infallible.  Approximately  15%  of  palpable  breast  cancers  are  not 
detectable  with  mammography,  and  this  number  is  likely  higher  in  younger  women.'’  A 
combination  of  clinical  palpation  with  either  mammography  or  sonography  has  been 
shown  to  significantly  increase  the  sensitivity  and  specificity  of  breast  cancer  detection.^ 
One  of  the  greatest  difficulties  in  mammography  is  imaging  the  radiographically  dense 
breast.  Unfortunately,  women  with  mammographically-dense  breasts  have  a  risk  of 
breast  cancer  that  is  1.8 — 6.0  times  greater  than  that  of  women  the  same  age  with  little  or 
no  mammographic  density.^  Small  lesions  become  much  more  difficult  to  detect  when 
obscured  by  dense  coimective  tissues  and  ducts.  Several  recent  studies  have  demonstrated 
that  sonography  has  higher  sensitivity  for  breast  cancer  detection  than  mammography 
alone,  ’  mammography  combined  with  physical  examination,’’*  or  mammoscintigraphy.® 
Furthermore,  the  sensitivi^  of  mammography  decreases  significantly  with  increasing 
mammographic  density.’’’”  Hormone  replacement  therapy  reduces  the  sensitivity  of  x-ray 
mammography’  ’  and  increases  the  need  for  alternate  diagnostic  tools. 

In  an  effort  to  improve  the  sensitivity  of  palpation  and  provide  quantitative 
measures  of  ‘palpable,’  research  groups  around  the  world  are  actively  working  toward 
imaging  technologies  that  display  quantitative  maps  of ‘tissue  stiffiiess.’  This  report 
describes  the  physics  of  palpation  and  uses  that  information  to  xmderstand  the  limitations 
of  palpation.  That  basic  physical  understanding  is  then  used  to  describe  the  various 
approaches  to  these  imaging  technologies.  The  emphasis  then  turns  to  elasticity  imaging 
systems  and  the  development  of  an  elasticity  imaging  system  that  is  implemented  on 
commercial  sonography  system  and  displays  real-time  elasticity  images  with  freehand 
scanning.  Results  of  preliminary  tests  of  the  utility  of  that  system  for  diagnosing  breast 
abnormalities  are  then  described. 

Previous  reviews  of  elasticity  imaging  using  ultrasound  are  available.’^’’^  This 
report  updates  those  prior  reviews  and  emphasizes  a  specific  real-time  elasticity  imaging 
system  and  results  obtained  with  that  system. 
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The  physics  of  palpation 

An  understanding  of  how  palpation  works  can  be  obtained  by  examining  the  basic 
physics  of  applying  an  external  deformation  to  an  object.  Begin  with  a  simple  model  of 
forces  and  deformation.  A  standard  concept  presented  in  introductory  physics  is  the 
elastic  deformation  of  a  simple  spring  (a  1-dimensional  object)  due  to  a  Imown  applied 
force.  Figure  1  illustrates  the  typical  simple  experiment  to  study  elasticity.  A  known  mass 
suspended  from  a  simple  spring  results  in  a  measurable  elongation  of  that  spring. 
Suspending  a  different  known  mass  results  in  a  different  elongation  of  the  spring.  Each 
mass  in  the  standard  gravitational  field  of  the  earth  places  a  known  force  on  the  spring. 
The  difference  in  these  forces  and  the  difference  in  the  elongations  of  the  spring  due  to 
those  forces  can  be  combined  using  Hooke’s  law  to  estimate  the  spring  constant,  k,  which 
is  characteristic  of  the  spring  and  quantifies  its  ‘stiffiiess.’ 

To  extend  the  concept  of  force  and  deformation  to  a  3-dimensional  (3-D)  object, 
consider  separately  the  forces  and  resulting  displacements.  The  analysis  will  be 
simplified  by  assuming  that  the  material  is  homogeneous  and  isotropic  (meaning  that  the 
material  properties  are  uniform  m  composition  without  any  directional  dependence  in 
elasticity).  Ignore  the  class  of  forces,  called  ‘body  forces,’  which  act  on  all  volume 
elements  of  the  material  (such  as  gravity  and  inertia).  The  class  of  forces  to  consider  is 
called  ‘surface  forces’  because  they  have  imits  of  force  per  xmit  area  and  can  be  viewed  as 
acting  on  a  surface  element  of  the  object.  That  surface  element  is  not  necessarily  on  the 
exterior  boundary  of  the  object,  but  can  be  a  surface  of  an  arbitrary  interior  volume 
element.  The  orientation  of  that  surface  is  described  by  a  vector  that  is  perpendicular  to 
the  surface  element  (a  normal  vector),  thus  a  3-D  coordinate  system  (x;,  i=l,2,3  or 
xi,X2,X3)  is  required  to  describe  the  normal  vector.  A  force  acting  on  that  surface  element 
has  a  magnitude  and  direction  (force  is  a  vector  quantity),  and  the  direction  of  that  force 
is  not  necessarily  perpendicular  to  the  surface  element.  Thus,  to  describe  the  direction  of 
the  force  vector  also  requires  a  3-dimensional  coordinate  system  (yj,  j=l,2,3).  To 
maintain  generality  and  simplicity  (to  obtain  principle  components)  in  the  description  of 
the  surface  force  two  separate  3-D  coordinate  systems  are  used  (xj  and  yj).  Collapsing  the 
arbitrary  surface  element  to  a  point  we  obtain  a  ‘stress  tensor.’  A  tensor  is  a 
generalization  of  the  concept  of  a  vector;  tensor  calculus  is  used  to  study  the  derivatives 
of  vector  fields.  The  stress  tensor,  ajj,  is  a  3x3  matrix  corresponding  to  the  nine 
combinations  available  by  combining  the  two  independent  3-D  coordinate  systems  of  the 
force  and  the  surface  element  on  which  it  acts. 

Similarly,  consider  the  displacement  of  a  volume  element  acted  on  by  an  external 
force.  If  the  motion  does  not  involve  a  change  of  volxxme  or  shape  of  the  object,  the 
motion  is  termed  ‘rigid  body  motion.’  If,  on  the  other  hand,  the  object  is  deformed 
(changes  shape  or  volume)  as  a  result  of  the  external  force  the  description  of  motion  is 
again  more  complex.  A  3-D  coordinate  system  is  required  to  describe  the  motion  in 
space.  To  maintain  generality  and  simplicity  (to  obtain  principle  components)  in  the 
description  of  the  gradients  (spatial  rate  of  change)  of  that  deformation  (strain  is  the 
spatial  rate  of  change  of  displacement)  another  3-D  coordinate  system  is  required.  The 
strain  tensor,  8^,  is  thus  another  3x3  matrix  corresponding  to  the  nine  combinations 
available  by  combining  these  two  independent  3-D  coordinate  systems. 

An  equation  to  relate  the  nine-component  stress  tensor  to  the  nine-component 
strain  tensor  is  called  a  constitutive  equation.  The  form  of  the  constitutive  equation 
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depends  on  whether  a  material  is  a  fluid  (an  ideal  fluid  with  no  viscosity  or  a  Newtonian 
viscous  fluid),  purely  elastic  (e.g.,  an  idealized  solid)  or  viscoelastic  (neither  purely 
viscous  nor  purely  elastic).  In  a  purely  elastic  (lossless)  deformation,  the  stress  is 
dependent  only  on  the  strain 

^ij  —  Cjjld  Sic]  • 

This  equation  is  analogous  to  Hooke’s  law  for  the  1-D  spring,  but  it  accoimts  for  forces 
and  deformations  in  all  three  directions.  The  quantity  Cyki  is  the  ‘modulus  tensor’  of 
elastic  coefficients  and  is  the  equivalent  of  the  spring  constant,  k,  used  to  describe  the 
deformation  of  a  spring.  The  four  subscripts  indicate  that  four  sets  of  3-D  coordinate 
systems  are  required  for  a  general  description  of  the  relationship  between  the  stress  and 
strain  tensors,  and  thus,  Cyki  has  (3"*)  81  components.  The  stress  and  strain  tensors  are 
symmetric  and  therefore  each  contains  at  most  six  independent  components.  Therefore, 
the  modulus  tensor  for  infinitesimal  elastic  deformations  is  also  symmetric  and  contains 
at  most  36  independent  components.  It  can  be  shown*'*  that  by  assiuning  a  material  to  be 
completely  isotropic,  the  number  of  independent  elastic  coefficients  is  reduced  to  two 
(called  the  Lame  constants).  A  more  detailed  description  of  stress  and  strain  can  be  found 
in  any  text  on  continuum  mechanics  (see,  for  example,  reference  15). 

The  elastic  coefficients  that  describe  the  behavior  of  a  material  are  absolute 
measures  of  intrinsic  properties  of  the  material.  Estimating  these  quantities  requires 
measurements  of  stresses  and  strains  under  well-characterized  experimental  conditions. 
For  example,  the  viscoelastic  properties  of  many  soft  tissues  under  cyclic  uniaxial 
loading  are  found  to  depend  on  the  strain  range,  the  strain  rate,  measurement  temperature, 
etc.*^  It  is  often  easier  to  simplify  the  experiment  and  only  measure  components  of  the 
surface  stress  distribution  or  the  internal  strain  distribution.  The  drawback  is  that  stress  or 
strain  alone  are  relative  quantities  and  are  not  intrinsic  to  the  material  under  study. 

The  basic  physics  of  elasticity  (stress  and  strain)  can  be  used  to  understand  the 
limitations  of  palpation.  Engineers  often  use  a  computational  tool  called  ‘finite  element 
analysis’  (FEA)  to  study  the  behavior  of  objects  under  external  forces  or  deformations. 
FEA  was  used  to  simulate  the  stress  and  strain  involved  when  deforming  a  uniform  block 
containing  a  spherical  inclusion,  as  shown  in  Figure  2.  The  iq)per  surface  of  the  block  is 
uniformly  displaced  by  1%  of  the  total  height.  The  lower  surface  is  allowed  to  move 
freely  laterally  and  the  sides  have  unrestricted  motion.  The  simulation  shows  the 
distribution  of  stress  and  profiles  of  that  stress  distribution  are  plotted  on  the  right.  In 
palpation,  the  fingers  press  on  the  tissue  to  deform  it  and  then  sense  the  stress  distribution 
that  results.  The  simulation  shows  that  as  the  sphere  moves  further  away  from  the  surface 
(profile  further  from  the  sphere),  the  variation  in  stress  across  the  profile  decreases 
suggesting  that  the  sphere  would  be  more  difficult  to  palpate  (less  stress  contrast 
available  to  the  fingers  to  sense)  as  it  is  placed  deeper  in  the  block. 

In  vitro  tissue  studies 

The  most  common  approach  to  studying  the  viscoelastic  properties  of  soft  tissues 
is  to  sinusoidally  deform  in  vitro  samples  of  tissue,  measure  the  force  required  to  induce 
the  deformation  and  study  the  phase  relationship  between  force  and  displacement.  In 
vitro  studies  of  the  viscoelastic  properties  of  breast  tissue*^’*’  have  demonstrated  several 
findings  that  are  significant  to  elasticity  imaging  (see  Table  1).  First,  for  cyclic  load- 
unload  experiments,  there  is  little  phase  delay  between  the  sinusoidal  deformation  and 
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response  (strain  and  stress)  for  compression  frequencies  near  IHz.  This  shows  that  the 
energy  required  to  deform  the  tissue  is  nearly  completely  recovered  when  the  deforming 
force  is  released  (nearly  lossless  deformation).  Thus,  in  vitro  breast  tissue  behaves  as  a 
nearly  completely  elastic  medium  at  these  strain  rates,  and  the  viscous  component  can  be 
ignored.  These  deformation  motion  frequencies  are  typical  of  that  used  in  clinical 
ultrasound  breast  exams  with  compression.  Second,  the  stress-strain  relationship  for  most 
breast  tissues  is  non-linear,  and  the  degree  of  non-linearity  varies  with  tissue  type. 
(Materials  with  linear  stress-strain  relationships  exhibit  stress  that  is  directly  proportional 
to  strain;  that  is,  they  exhibit  constant  ‘stiffiiess.’  Materials  with  nonlinear  stress-strain 
relationships  change  stif&iess,  most  commonly  getting  stiffer,  as  they  are  deformed.) 
Third,  the  elastic  moduli  of  breast  tissue,  obtained  from  the  slope  of  the  stress-strain 
curves,  vary  significantly  among  breast  tissue  types  and  strain  range.  In  summary,  breast 
tissue  is  mostly  elastic  for  the  strain  rates  likely  encountered  with  freehand  scanning, 
object  contrast  is  likely  high  in  strain  and  modulus  images,  contrast  will  likely  be 
different  for  different  lesion  types,  and  contrast  will  likely  change  with  increasing 
compression. 


Imaging  the  elastic  properties  of  tissue 

Approaches  to  elasticity  imaging  can  be  classified  by  the  modality  of  the  signal 
source  (primarily  ultrasound  or  MRI),  the  mechanical  parameter  estimated  (e.g.,  stress, 
strain,  or  modulus),  or  a  descriptor  of  the  experimental  procedme  (“dynamic”  or  “(quasi- 
)static”  techniques).  The  mechanical  properties  estimated  with  these  techniques  are 
related.  As  described  above,  stress  and  strain  are  mutually  responsive  quantities,  but  they 
are  not  intrinsic  material  properties.  Images  of  stress  and  strain  are  maps  of  a  parameter 
relative  to  its  surroundings  (as  a  mammogram  maps  the  relative  x-ray  attenuation,  for 
example).  Elastic  moduli  are  intrinsic  material  properties  generally  described  with  a 
matrix  (as  described  above),  but  for  practicality  experimental  conditions  are  manipulated 
and  material  properties  (such  as  incompressibility,  homogeneity  and  isotropy)  are 
assumed  so  that  the  size  of  this  matrix  is  reduced  to  one  or  two  parameters. 

Several  research  groups  are  developing  techniques  for  imaging  the  stress 
distribution.  Most  notable  among  the  stress  imaging  techniques  is  the  work  from 
Wellman,  et  al.,**  using  a  piezo-resistive  sensor  array  (TekScan,  Inc.,  Boston,  MA) 
coupled  to  a  position  tracking  system.  This  system  closely  mimics  the  mechanics  of 
palpation  and  demonstrates  a  strong  correlation  between  the  size  of  the  lesion  measured 
with  the  tactile  system  and  the  lesion  size  measured  following  resection.  The 
performance  for  small  lesions  (less  than  10mm)  that  are  relatively  deep  (more  than 
10mm)  remains  to  be  seen.  Also  noteworthy  is  the  work  of  Sarvazyan*^  in  which  he 
attempts  to  solve  the  inverse  problem  of  determining  the  3-D  modulus  distribution  that 
causes  the  measured  surface  pressure  distribution. 

Strain  imaging  has  received  the  most  attention  in  elasticity  imaging.  The  earliest 
implementations  used  M-mode  acquisition  and  cross  correlation  to  track  tissue  motion 
and  study  tissue  elasticity. In  later  studies  Doppler  processing  techniques  were  used 
to  track  differences  in  motion^^’^^  and  “sonoelasticity  imaging”  soon  followed.^'*  The 
Doppler  processing  techniques  were  the  first  “dynamic  techniques”  and  derived  their  data 
from  ultr^lsound.  ‘Static  compression  elastography’  is  the  most  common  approach  to 
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strain  imaging.  Numerous  groups  are  pursuing  ultrasound-based  strain  imaging  with 
efforts  in  algorithm  development,  ’  >  ■  ■  performance  evaluation,  ’  ’  and  clinical 
testing^^’^''  (representative  citations). 

The  basic  information  derived  in  strain  imaging  techniques  is  the  relative  tissue 
displacement.  An  imaging  system  (typically  ultrasound  or  MRI)  acquires  (pre¬ 
deformation)  data  corresponding  to  a  map  of  tissue  anatomy.  A  small  deformation  is 
applied,  either  through  an  external  compressor  or  physiological  function  (breathing, 
cardiac  pressure  variations,  etc.),  and  another  (post-deformation)  map  of  the  anatomy  is 
acquired.  The  displacement  field  in  the  deformed  tissue  is  estimated  by  comparing  these 
two  maps  of  anatomy.  Mechanical  strain  is  estimated  by  calculating  the  gradient  (the 
spatial  rate  of  change)  of  the  displacement  field.  In  ultrasound  the  displacement  along  the 
acoustic  beam  propagation  (axial)  direction  can  be  estimated  far  more  accxurately  and 
with  higher  precision  than  in  the  lateral  or  elevational  directions.^^ 

An  important  aspect  for  clinical  acceptance  of  ultrasound  strain  imaging  is  the 
technique  for  deforming  the  soft  tissue  between  image  pairs.  Most  phantom  experiments 
in  the  literature  used  motorized  compression  devices  and  extensive  fixtures,  liiese 
devices  are  not  likely  to  gain  clinical  acceptance  because  they  either  limit  the  locations 
that  strain  imaging  can  be  applied  or  are  time  consuming  to  incorporate.  Freehand 
scanning,  where  tissue  is  deformed  with  the  surface  of  the  transducer,  is  desirable.^^’^'*’^* 

Developing  a  real-time  strain  imaging  system  that  allows  freehand  scanning  is 
essential  for  clinical  utilization.  The  strain  imaging  algorithm  must  be  computationally 
efficient,  insensitive  to  motion  irregularities  and  track  tissue  motion  in  2-D  (eventually  3- 
and  4-D).  Block  matching  (template  matching)  algorithms  are  widely  used  in  image 
processing  applications  for  tracking  motion.  The  most  notable  application  is  movie  image 
compression  algorithms  such  as  MPEG.  The  use  of  block  matching  in  ultrasonic  imaging 
was  first  reported  by  Trahey,  et  al.  for  blood  flow  estimation.^®  Block  matching  is  a  good 
candidate  since  it  is  simple  in  principle  and  is  capable  of  tracking  motion  in  2-D. 
However,  for  strain  imaging,  the  algorithm  needs  to  be  modified  to  increase  its 
computational  efficiency  and  insensitivity  to  decorrelation  noise.^^  (Decorrelation  is  a 
measure  of  how  similar  two  signals  are.  That  similarity  is  measured  with  cross  correlation 
or  surrogate  measures  of  correlation.  Echo  signals  decorrelate  when  there  is  high 
electronic  noise  or  when  there  is  large  deformation  of  the  tissue.) 

There  has  been  less  attention  focused  toward  strain  imaging  systems  than  toward 
strain  imaging  algorithms,  data  simulation,  and  performance  testing.  Bamber,  et  al.,  have 
reported^®  their  progress  in  freehand  elasticity  imaging.  Their  system  lacked  real-time 
feedback  in  the  data  acquisition  process,  but  still,  tiiey  found  that  it  is  possible  to  obtain 
good  elasticity  data  with  freehand  scanning.  Their  rate  of  success  was  relatively  low,  and 
significant  pre-  and  post-processing  was  necessary  to  obtain  accurate  displacement 
.  estimates.  The  system  reported  by  Garra,  et  al.^^  employed  a  modified  mammography 
paddle  with  a  hole  cut  out  to  provide  an  acoustic  window.  This  allowed  (relatively)  easy 
correlation  with  the  mammogram.  However,  the  acoustic  data  acquisition  system  was 
crude.  The  system  only  allowed  scanning  with  a  5MHz  transducer — lower  than  the 
standard  of  its  day  (7.5MHz) — ^and  had  significantly  poorer  performance  than  current 
systems.  In  addition,  the  digitization  was  external  to  the  ultrasound  scanner  resulting  in 
reduced  electronic  SNR  and  increased  timing  jitter  in  the  acquired  echo  signals.  The 
increased  jitter  significantly  reduces  the  performance  of  displacement  estimates  in  strain 
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imaging.  The  current  system  used  by  that  group  incorporates  a  mid-range  ultrasound 
scanner  with  a  5-axis  motor  controlled  compression  system.^^  The  first  real-time 
elasticity  imaging  system  was  developed  for  prostate  imaging.^*  Data  were  acquired  in  a 
sector-shaped  scan  from  an  endo-cavity  transducer  and  a  1-D  tracking  method  was  used. 
As  a  result,  elasticity  image  frame  rates  were  quite  high  at  the  expense  of  image  quality. 

The  in  vivo  studies  of  strain  imaging  reported  by  Garra,  et  al.,^^  demonstrated 
that  strain  imaging  has  merit  in  differentiating  among  solid  tumors  in  breasts.  Their  most 
significant  finding  was  that  invasive  ductal  carcinomas  are  significantly  wider  in  strain 
images  than  in  the  corresponding  B-mode  image,  and  this  difference  is  likely  due  to  the 
desmoplastic  reaction  that  surroimds  this  tumor  type. 

Modulus  imaging  has  also  been  investigated,  and  there  are  three  primary 
approaches  in  the  literature.  The  first  approach  estimates  the  shear  wavelength  in  tissue 
and  from  this  directly  estimates  the  shear  modulus  of  the  tissue.^®’'”’’'"  The  other 
techniques  require  simultaneous  measurements  of  stress  and  strain  and  require 
assumptions  regarding  the  bormdaiy  conditions  of  the  experiment.'*^’^^’'*^’^  ’'’^Compared  to 
strain  imaging,  modulus  imaging  has  lower  spatial  resolution,  higher  noise,  and  the 
assumptions  regarding  boundary  conditions  can  result  in  biased  estimates.  However, 
estimating  an  intrinsic  tissue  parameter,  instead  of  the  relative  parameters  estimated  in 
stress  or  strain  images,  makes  this  an  attractive  approach. 

There  are  also  methods  under  development  that  use  acoustic  radiation  force  to 
deform  tissue  and  study  tissue  viscoelasticity^’’'*^’'*^  with  promising  results.  Other  novel 
approaches  to  describing  the  viscoelastic  behavior  of  tissues,  such  as  those  reported  by 
Greenleaf,  et  al.,^®’^'  are  also  under  investigation. 

Early  work  in  strain  imaging  demonstrated  the  limitations  of  tracking  motion  in  1- 
D  and  motivated  the  development  of  2-D  and  3-D  motion  tracking  algorithms  for 
elasticity  imaging.  ’  ’  Those  studies  demonstrated  that  1-D  tracking  failed  to  correctly 
track  motion  in  a  3cm  wide  field  of  view  with  as  little  as  0.6%  compression,  and  motion 
tracking  errors  got  increasingly  worse  with  increased  compression.  However,  by  using  2- 
D  tracking  algorithms  that  appropriately  compensate  for  lateral  motion,  high  contrast-to- 
noise  images  of  mechanical  strain  could  be  obtained  with  compressions  of  more  than  5% 
in  phantoms  (see  Figure  3).  The  basic  approach,  called  ‘companding,’  was  to  use  2-D 
motion  tracking  to  align  (warp)  either  the  pre-  or  post-deformation  data  field  prior  to  1-D 
cross  correlation. 

Other  early  work  also  demonstrated  the  need  to  control  motion  during  elasticity 
imaging  experiments.  The  images  in  Figure  4  demonstrate  that  it  is  essential  to  control 
the  motion  during  deformation  especially  with  regard  to  elevation  motion.  A  typical 
clinical  ultrasound  imaging  system  acquires  echo  data,  nominally,  from  a  plane  of  tissue. 
Any  out-of-plane  motion  of  tissue  will  result  in  echo  signal  decorrelation  and  reduced 
elasticity  image  quality. 

The  key  to  obtaining  high  quality  elasticity  images  is  the  quality  of  the  motion 
tracking  algorithm.  Ultrasound  radio-frequency  (rf)  echo  signals,  the  same  data  used  to 
form  a  B-mode  image,  are  used  as  a  map  of  anatomy.  Those  same  signals  are  used  to 
track  the  deformation  of  the  anatomy.  The  task  is  to  accurately  track  the  anatomical 
deformation  with  minimal  uncertainty  (displacement  estimate  variance  or  covariance).  A 
review  of  many  of  the  techniques  used  for  tracking  tissue  motion  with  ultrasound  can  be 
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found  in  reference  53.  A  tutorial  on  the  general  topic  of  waveform  coherence  and  time- 
delay  estimation  can  be  found  in  reference  54. 

A  review  the  assumptions  used  in  signal  correlation  analysis  can  help  to 
appreciate  the  difference  between  many  motion  tracking  algorithms.  A  typical 
assumption  in  motion  tracking  based  on  time  delay  of  ultrasound  echo  signals  is  that  the 
deformation  of  the  tissue  is  minimal  (or  recoverable)  within  the  echo  signal  segment 
being  tracked.  Another  common  assumption  is  that  the  observation  window  (data 
segment  length)  is  large  compared  to  the  time  delay.  Thus,  a  relatively  long  data  segment 
is  needed  to  avoid  ambiguous  displacement  estimates  (referred  to  as  ‘peak  hopping’).  The 
plot  in  Figure  5  demonstrates  that  with  an  if  echo  segment  as  short  as  3mm  and  with  only 
1.5%  axial  strain,  there  is  obvious  echo  signal  decorrelation  between  the  pre-  and  post¬ 
deformation  A-lines.  However,  the  single  large  peak  in  the  cross  correlation  function 
plotted  in  Figure  5  demonstrates  that  there  is  little  ambiguity  in  the  time  delay  required  to 
match  pre-and  post-deformation  signals.  Using  shorter  rf  echo  segments  in  motion 
tracking  reduces  the  decorrelation  within  the  echo  signal  segment,  and  increases  the 
waveform  coherence,  as  shown  in  Figure  6.  However,  short  data  segments  increase  the 
likelihood  of  time  delay  ambiguity  (e.g.,  a  one-wavelength  segment  of  rf  looks  very  much 
like  many  other  one-wavelength  segments).  Using  multiple  (usually  adjacent)  A-line 
segments  reduces  the  likelihood  of  ambiguity,  as  shown  in  Figure  7.  Short  data  segments 
also  demonstrate  the  benefit  of  interpolating  (up-sampling)  the  rf  echo  signal.  The 
waveforms  shown  in  Figure  6  illustrate  that  waveform  coherence  would  improve  if  time- 
delays  of  less  than  one  sample  were  available.  An  alternative  is  to  interpolate  the 
correlation  function,  but  this  requires  a  model  for  the  functional  form  of  the  cross 
correlation  function.  If  up-sampling  the  rf  echo  signal  can  be  justified,  it  reduces  the  need 
for  an  accurate  model  of  the  cross  correlation  function  when  interpolating  sub-sample 
displacement  estimates,  as  show  in  Figure  8. 

Development  of  a  real-time  strain  imaging  system 

Experience  in  developing  motion  tracking  algorithms  and  experiments  with 
phantoms  and  in  vitro  tissues  suggest  criteria  for  a  clinically  viable  elasticity  imaging 
system.  First,  the  system  must  track  tissue  motion  in  2-D  (or  3-D,  if  available)  for  high 
contrast-to-noise  images.^’’^*  Second,  the  system  should  use  short  2-D  data  segments 
(kernels)  for  motion  tracking  to  minimize  decorrelation  within  the  data  segments  and  to 
minimize  time-delay  ambiguity.  Third,  the  system  should  provide  real-time  elasticity 
images,  as  well  as  normal  B-mode  images,  to  allow  the  user  to  monitor  the  images  being 
acquired  and  manipulate  the  transducer  array  with  fi:eehand  scanning  thus  insuring  that 
the  tissue  motion  is  suitable  for  forming  high-quality  elasticity  images.  In  addition,  the 
data  acquisition  technique  should  be  similar  to  that  currently  used  in  sonography  to 
increase  the  likelihood  of  clinical  acceptance.  A  large  deviation  firom  standard  clinical 
practice  would  likely  receive  a  more  skeptical  assessment  by  potential  users  than  a  subtle 
modification  to  current  practice. 

A  novel  motion  tracking  algorithm  has  been  developed  and  implemented  on  a 
clinical  ultrasound  imaging  system  (SONOLINE  Elegra,  Siemens  Medical  Solutions, 
Issaquah,  WA).^^  Phase-sensitive  (I-Q)  echo  data  are  processed  internally  in  real-time  on 
the  Elegra  to  estimate  displacement  and  strain.  The  system  can  use  any  of  the  linear  array 
transducers  available  on  the  Elegra  and  is  compatible  with  Tissue  Harmonic  Imaging  on 


10 


that  system.  The  system  displays  B-mode  and  strain  images  side-by-side  on  the  normal 
system  display  at  about  7  frames/sec.  A  region  of  interest  (ROI)  is  displayed  in  the  B- 
mode  image  and  displacement  and  strain  are  estimated  for  tissue  within  that  ROI.  The 
size  and  location  of  the  ROI  can  be  manipulated  with  front  panel  controls.  When 
scaiming,  the  normal  freeze  and  cine  capabilities  of  the  system  are  available.  When  a 
sequence  of  data  is  acquired  and  stored  (frozen),  on-line  post-processing  capabilities 
allow  the  ROI  location  and  size  to  be  modified,  and  other  common  tools  such  as 
modifying  the  grayscale  mapping  are  available.  Initial  tests  of  the  elasticity  image  noise 
and  spatial  resolution  are  found  in  reference  29.  Spherical  lesions  as  small  as  2.4nun 
diameter  that  are  three  times  stiffer  than  the  background  were  easily  displayed.  The 
protocol  for  clinical  testing  of  this  system  was  approved  by  the  Humans  Subjects 
Committee  at  the  University  of  Kansas  Medical  Center  where  that  initial  work  was 
performed. 

A  critical  issue  in  the  development  and  utility  of  any  imaging  system  is  the 
achievable  spatial  resolution  for  a  given  task.  The  ability  to  image  a  3mm  diameter 
sphere  in  a  phantom  is  encouraging.  More  importantly,  those  encouraging  results  are 
corroborated  with  the  ability  to  image  small  structures  in  vivo.  For  example,  images  of  an 
in  vivo  3mm  cyst  are  shown  in  Figure  9.  Although  the  ability  to  image  small  structures  in 
vivo  is  clearly  demonstrated,  the  required  contrast  to  view  objects  of  a  specific  size  is 
unknown.  Investigations  are  underway  to  evaluate  this  through  contrast-detail  analysis.^^ 

The  ability  to  acquire  and  view  long  sequences  of  elasticity  images  has  provided 
the  opportunity  to  observe  nonlinear  elastic  behavior  of  in  vivo  tissues.  Nonlinearity  in 
the  stress-strain  relationship  of  tissue  was  observed  with  in  vitro  breast  tissues,*’  and  was 
therefore  expected  with  in  vivo  tissues,  but  was  only  recently  observed  with  the 
availability  of  a  real-time  elasticity  imaging  system.^'*  Figure  10  shows  an  example  of  the 
implications  of  nonlinear  elasticity  in  strain  imaging.  At  low  preload  (transducer  barely  in 
contact  with  the  skin  surface  and  applying  minimal  pressure)  the  fibroadenoma  appears 
dark  in  the  strain  image.  As  the  preload  is  increased  (pressure  applied  with  the  transducer 
increasing  deformation),  the  strain  image  contrast  of  the  fibroadenoma  (its  stiffness 
relative  to  the  surroimding  tissue)  decreases.  This  behavior  might  explain  why  others 
have  found  that  some  fibroadenomas  are  not  visible  in  single  strain  images.^^ 

One  of  the  significant  findings  in  prior  clinical  trials  of  in  vivo  elasticity 
imaging^^  was  that  the  size  of  a  breast  lesion  displayed  in  strain  images,  relative  to  it  size 
in  a  normal  B-mode  image,  appears  to  be  a  significant  criterion  for  differentiating 
malignant  from  benign  breast  lesions.  Figures  1 1  and  12  show  examples  of  the  B-mode 
and  strain  image  pair  for  a  fibroadenoma  and  an  invasive  ductal  carcinoma,  respectively. 
In  each  case,  the  lesion  is  traced  in  the  B-mode  image  and  that  tracing  is  reproduced  in 
the  respective  strain  image.  The  lesion  boimdary  traced  for  benign  lesions  have  about  the 
same  size  and  shape  in  the  two  image  types.  However,  the  lesion  boundary  traced  in  B- 
mode  images  of  invasive  ductal  carcinomas  is  much  smaller  than  the  lesiori  displayed  in 
the  respective  strain  image.  On  average,  the  area  of  these  carcinomas  displayed  in  strain 
images  is  three  times  larger  than  that  of  B-mode  images.^"*  It  is  postulated^^’^  that  the 
increased  size  of  carciuomas  in  strain  images  is  due  to  the  desmoplasia  that  often 
sxuTOunds  invasive  ductal  carcinoma. 

To  test  the  utility  of  relative  lesion  size  for  differentiating  between  benign  and 
malignant  breast  lesions,  five  observers  individually  viewed  a  set  of  image  sequences 
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from  in  vivo  breast  elasticity  imaging.  Each  observer  selected  the  image  pair  from  a 
sequence  (movie  loop)  that  was  most  representative  of  the  B-mode  and  strain  image  pair 
from  that  sequence.  Each  observer  then  traced  the  outline  of  the  lesion  in  each  image  type 
and  measured  the  width  and  height  of  the  lesion  in  each  image.  This  was  repeated  for 
data  from  97  movie  loops  of  55  unique  lesions  from  29  patients.  A  plot  of  the  average 
lesion  area  for  each  lesion  measured  by  the  group  of  observers  is  shown  in  Figure  13. 
These  data  are  consistent  with  those  reported  by  others^^  and  suggest  that  elasticity 
imaging  may  be  a  useful  tool  to  improve  the  utility  of  breast  sonography.  If  the  ratio  of 
lesion  size  in  strain  images  versus  B-mode  images  proves  to  be  a  sensitive  criterion  for 
increasing  confidence  of  a  benign  diagnosis,  the  fraction  of  biopsies  that  prove  to  be 
benign  tissue  will  likely  be  reduced  at  significant  savings  in  health  care  expense  and 
trauma  to  the  patient  and  their  family  and  friends. 

Summary 

Elasticity  imaging  is  a  relatively  new  technique  for  studying  the  stifi&iess  of 
tissue.  The  information  acquired  with  these  techniques  is  similar  to  that  obtained  with 
manual  palpation,  but  elasticity  imaging  is  more  sensitive  and  less  subjective  than 
palpation.  Further,  the  information  is  provided  in  an  image  format  so  that  it  can  be 
compared  with  data  from  other  image  modalities  and  it  can  more  easily  be  documented 
and  shared  with  others. 

Several  interesting  approaches  to  elasticity  imaging  are  currently  being 
investigated  by  research  groups  around  the  world.  Different  approaches  provide  different 
information  about  the  viscoelastic  properties  of  tissue.  Many  of  these  approaches 
emphasize  the  elastic  properties  of  tissue  due  to  the  techniques  of  data  acquisition. 

At  least  one  method  for  elasticity  imaging  is  under  development  that  produces 
images  of  mechanical  strain  in  real-time  using  a  freehand  scanning  technique  very  similar 
to  that  of  standard  breast  sonography  examinations.  The  system  is  integrated  into  a 
clinical  sonography  system  without  any  external  equipment  and  involves  software 
changes  only.  In  vivo  tests  of  this  system  have  demonstrated  the  ability  to  image  small 
breast  lesions  with  confidence.  It  has  also  allowed  the  visualization  the  effects  of 
nonlinear  elasticity  of  in  vivo  breast  tissues.  Further  investigations  with  this  system 
suggest  that  benign  breast  lesions  are  about  the  same  size  and  shape  in  B-mode  and  strain 
images,  but  invasive  ductal  carcinomas  tend  to  be  significantly  larger  in  strain  images 
than  in  the  corresponding  B-mode  images.  This  fact  suggests  that  elasticity  imaging 
might  increase  the  utility  of  breast  sonography  and  might  be  offered  in  clinical 
sonography  systems  in  the  near  future. 
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Lesion  type 

5%  Precompression 
10%/s  Strain  Rate 

20%  Precompression 
20%/s  Strain  Rate 

Normal  Fat 

19±7 

20±6 

Normal  Glandular 

33  ±11 

57  ±19 

Fibrous 

107±32 

233  ±  59 

Ductal  Carcinoma  In  Situ 

25  ±4 

301  ±58 

Invasive  Ductal  Carcinoma 

93  ±33 

490  ±112 

Table  1.  Results  of  measuring  the  elastic  modulus  of  in  vitro  breast  tissue  samples  for  a 
variety  of  breast  tissue  types  with  two  different  strain  ranges  and  strain  rates.  Ductal 
carcinoma  has  the  largest  difference  for  modulus  measurements  in  these  two  conditions 
suggesting  it  has  the  most  nonlinear  stress-strain  relationship  among  tissues  measured. 
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F  =  m  a 


AF  =  k  Ax 

Spring  Constant  k 


Figure  1.  Introductory  mechanics  describes  the  behavior  of  a  spring  supporting  different 
masses.  The  diagram  illustrates  how  one  spring  would  be  stretched  to  two  different 
lengths  by  two  different  masses.  Hooke’s  law  describes  this  behavior  and  can  be  used  to 
characterize  the  spring  under  small  deformations.  Newton’s  second  law  equates  a  force, 
F,  with  a  mass,  m,  and  an  acceleration,  a.  A  known  mass  suspended  from  a  spring  exerts 
a  known  force  due  to  the  acceleration  of  gravity.  Hooke’s  law  relates  die  difference  in 
stretch  of  the  spring.  Ax,  due  to  the  change  of  force,  AF,  resulting  from  suspending 
different  masses.  The  proportionality  constant,  k,  characterizes  the  stiffiiess  of  the  spring. 
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Figure  2.  Palpation  can  be  approximated  with  a  simulation  tool  (FEA).  The  center  image 
is  the  axial  stress  distribution  resulting  from  a  uniform  displacement  of  the  top  surface  of 
a  block  containing  a  spherical  object.  The  plot  on  the  right  shows  profiles  of  &e  stress 
distribution  across  the  lines  in  the  center  image.  A  large  variation  in  the  stress  profile  is 
found  for  the  profile  close  to  the  spherical  object.  This  simulates  a  sphere  being  near  the 
surface  of  the  block.  As  the  profile  moves  further  away  from  the  sphere,  the  variation  in 
stress  across  the  profile  decreases  suggesting  the  sphere  would  be  more  difficult  to 
palpate  as  it  is  placed  deeper  in  the  block. 
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Figure  3.  Images  of  the  mechanical  strain  in  the  axial  direction  for  a  gelatin  phantom  with 
three  cylinders  which  are  three  times  stiffer  than  the  background.  The  top  row  shows 
strain  images  obtained  with  1-D  tracking  using  cross  correlation.  The  bottom  row  shows 
images  obtained  with  2-D  companding  (2-D  motion  tracking).  The  applied  deformation 
(top  and  bottom  rows),  from  left  to  right,  is  0.6%,  1.2%,  2.4%  3.6%,  4.8%  and  6.0% 
strain.  Lateral  expansion  (bulging)  occurring  with  axial  compression  causes  the  echo  A- 
lines  to  not  match  and  the  echo  signals  to  decorrelate  with  1-D  tracking.  2-D  tracking  is 
able  to  track  the  lateral  as  well  as  axial  motion  to  obtain  higher  contrast-to-noise  strain 
images. 
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Figure  4.  Images  acquired  near  the  edge  (in  elevation,  perpendicular  to  the  image  plane) 
of  a  gelatin  phantom  with  varying  deformation  (from  left  to  right,  0.6%,  1.2%  and  2.4% 
axial  strain).  The  phantom  was  bound  in  elevation  at  the  bottom  and  slid  freely  at  the  top. 
As  the  phantom  was  deformed,  the  top  slid  out  of  the  image  plane  in  elevation  resulting 
in  echo  signal  decorrelation.  Increasing  the  deformation  caused  greater  decorrelation. 
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Figure  5.  On  top  is  a  plot  of  the  pre-deformation  (red  dashed  line)  and  post-deformation 
(1.5%  axial  strain,  solid  blue  line)  for  a  3nun  segment  (140  rf  samples)  of  the  echo 
signals  from  the  center  of  a  gelatin  phantom.  The  post-deformation  signal  has  been 
shifted  in  time  to  match  the  pre-deformation  signal  as  well  as  possible.  The  deformation 
has  caused  decorrelation  in  the  echo  signal  which  reduces  coherence  (cross  correlation 
coefficient  =  0.87).  Below  is  a  plot  of  the  cross  correlation  function  comparing  the  pre- 
and  time-delayed  post-deformation  rf  echo  signals  shown  above.  The  single  large  positive 
peak  suggests  there  is  little  ambiguity  in  the  correct  delay  required  to  match  these  signals. 
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Figure  6.  On  top  is  a  plot  of  the  pre-deformation  (red  dashed  line)  and  post-deformation 
(1 .5%  axial  strain,  solid  blue  line)  for  a  0.24mm  segment  (1 1  rf  samples)  of  the  echo 
signals  from  the  center  of  a  gelatin  phantom.  The  post-deformation  signal  has  been 
shifted  in  time  to  match  the  pre-deformation  signal  as  well  as  possible.  Little 
decorrelation  in  the  echo  signals  within  this  short  echo  segment  results  in  high  coherence 
(cross  correlation  coefficient  =  0.96).  Below  is  a  plot  of  the  cross  correlation  function 
comparing  the  pre-  and  time-delayed  post-deformation  rf  echo  signals  shown  above. 
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Figure  7.  Plots  of  the  pre-deformation  (solid  blue)  and  the  time-shifted  post-deformation 
(1.5%  strain,  dashed  red)  rf  echo  signal  from  five  adjacent  A-lines  obtained  near  the 
center  of  a  gelatin  phantom.  Using  multiple  short  line  segments  reduces  decorrelation 
within  the  data  segment  while  simultaneously  reducing  the  ambiguity  in  time-delay 
estimation. 
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Figure  8.  In  the  left  column  are  plots  of  the  pre-deformation  (red  dashed  line)  and  post¬ 
deformation  (1.5%  axial  strain,  solid  blue  line)  for  a  0.24mm  segment  of  the  echo  signals 
acquired  at  36MHz  sampling  frequency  up-sampled  to  72MHz  (top  row)  and  144MHz 
(bottom  row).  The  post-deformation  signal  has  been  shifted  in  time  to  match  the  pre¬ 
deformation  signal  as  well  as  possible.  As  the  effective  sampling  interval  is  reduced,  the 
integer  time  delay  error  is  also  reduced  allowing  for  greater  waveform  coherence 
(correlation  coefficients  of 0.965  (Figure  6),  0.978,  and  0.998  were  obtained  with 
36MHz,  72MHz  and  144MHz  sampling,  respectively). 
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Figure  9.  B-mode  and  elasticity  (mechanical  strain)  images  of  a  3mm  diameter  in  vivo 
breast  cyst  demonstrating  that  small  in  vivo  structures  are  resolvable  in  strain  images. 
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Figure  10.  B-mode  and  strain  images  of  a  typical  fibroadenoma  at  two  different  amounts 
of  deformation  (preload).  For  both  elasticity  images  the  average  strain  in  the  image  is 
about  1.2%.  The  image  pair  on  the  left  was  acquired  with  the  ultrasound  transducer  just 
barely  in  contact  with  the  skin  surface.  At  low  preload  fibroadenomas  typically  are  stiff 
compared  to  their  surroimding  glandular  tissue.  As  the  preload  increases  the  strain  image 
contrast  of  the  fibroadenoma  decreases  and  the  lesion  becomes  nearly  equal  stifl&iess 
compared  to  the  surroimding  tissue. 
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Figure  1 1.  A  B-mode  and  strain  image  pair  of  a  fibroadenoma.  The  lesion  is  traced  in  the 
B-mode  image  and  that  tracing  is  displayed  in  the  strain  image.  The  lesion  size  and  shape 
in  the  two  image  types  are  very  similar. 
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Figure  12.  A  B-mode  and  strain  image  pair  of  a  scirrhous  invasive  ductal  carcinoma.  The 
lesion  is  traced  in  the  B-mode  image  and  that  tracing  is  displayed  in  the  strain  image.  The 
lesion  size  is  larger  in  the  strain  image  than  in  the  B-mode  image. 
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Figure  13.  A  plot  of  the  lesion  area  measured  in  strain  images  versus  the  same  lesion 
measured  in  the  corresponding  B-mode  images  for  cysts,  fibroadenomas  and  invasive 
ductal  carcinomas.  The  average  result  of  five  observers  is  plotted  and  error  bars  represent 
standard  deviations  among  those  measurements.  The  dashed  line  represents  equal  area  in 
both  image  types. 
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